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This  report  describes  the  design  of  the  Intelligent  Event  Identification  System  (ISEIS)  and  the  key  < 

processing  algorithms  of  the  system.  ISEIS  will  characterize  seismic  events  utilizing  discriminants  j 

(  to  classify  source  types,  explosion  or  earthquake,  and  by  comparing  waveform  features  with  j 

,  reference  events  to  estimate  how  similar  the  events  are  to  previously  observed  events.  Various  i 

graphics  interfaces  will  provide  explanations  to  the  analyst  about  the  event  characterization  results  | 

as  well  as  various  graphics  displays  of  the  data  processing  results.  ISEIS  will  initially  be  devel-  i 

oped  to  identify  small  events  at  regional  distances,  recorded  at  the  regional  arrays.  A  discrimination  j 

study  of  regional  discriminants  at  NORESS  for  mine  explosions  and  earthquakes  in  Norway  is  ! 

presented.  This  study  has  revealed  that  ratios  of  regional-phase  amplitudes  at  frequencies  above  6  ; 

Hz  provide  the  best  discrimination  capability.  Spectral  discriminants,  such  as  Lg  spectral  ratio,  ap-  1 

pear  to  be  less  useful.  However,  the  identification  of  time-independent  spectral  modulations  effec-  j 

tively  identifies  most  mine  explosions  which  are  ripple  fired.  A  process  called  Multiple  Event 
Recognition  System  (MERSY)  is  described  which  automatically  identifies  ripple-fired  explosions 
using  cepstral  methods.  Finally,  the  Dynamic  Time  Warp  (DTW)  algorithm  is  described  which  will  ! 
be  used  in  ISEIS  to  do  waveform  matching  of  regional  seismic  codas. 
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SUMMARY 


The  objective  of  this  project  is  to  design  and  develop  an  Intelligent  Event  Identification 
System,  or  ISEIS,  which  will  be  a  prototype  for  routine  event  identification  and  to  serve  as  a 
research  tool  for  discrimination  research.  In  stating  this  goal,  it  is  important  to  define  what  is 
meant  by  "intelligent."  We  have  tried  to  design  ISEIS  such  that,  when  presented  a  specific  seismic 
event,  the  system  will  approach  the  characterization  of  the  event  in  a  systematic  fashion  in  the  same 
way  that  a  seismologist  would.  Eor  example,  a  seismologist  would  be  aware  of  the  latest  results  of 
discrimination  research  and  would  try  to  apply  the  same  discriminants.  However,  the  seismologist 
would  be  careful  not  to  blindly  apply  the  discriminants,  since  they  may  be  affected  by  regional  path 
effects.  Thus,  it  is  important  that  event  identifications  be  "qualified"  based  on  what  the  limitations 
are  about  our  understanding  of  the  discriminants  in  terms  of  how  they  may  be  affected  by  non¬ 
source  effects.  The  analyst  should  also  be  flexible  enough  to  learn  new  discrimination  methods 
based  on  the  results  of  the  analysis  of  additional  data.  Finally,  any  seismologist  would  want  to 
compare  the  event  with  other  nearby  events  and/or  waveform  features  extracted  from  reference- 
event  waveforms.  Such  comparisons  would  be  important  in  any  informed  characterization  of  the 
event,  particularly  since  regional  seismic  discrimination  is  not  well  understood. 

ISEIS  will  be  capable  of  doing  this  kind  of  detailed,  systematic  analysis  cf  every  event 
which  is  detected  and  located  by  the  NMRD  system,  in  both  an  automated  and  interactive  mode.  In 
the  automated  mode,  the  system  will  apply  a  set  of  standard  discriminants  in  routine  fashion  on  any 
new'  events  which  appear  in  the  database.  In  the  interactive  mode,  the  user  can  review  the  results 
of  the  automated  processing.  Additional  analysis  can  be  done,  if  necessary,  to  improve  or  enhance 
the  processing  results.  Discriminants  can  then  be  run  again,  either  individually  or  all  together. 

This  project  has  been  a  combined  system  design  and  development  project  and  a  research 
project.  The  initial  ISEIS  system  will  be  an  "experimental  prototype,"  meaning  that  it  may  be 
modified  and  improved  in  the  future  after  it  has  been  evaluated.  As  of  this  writing,  the  first 
prototype  of  the  ISEIS  system  has  been  designed  and  development  is  nearing  completion. 


VII 


However,  the  system  itself  will  not  be  described  in  this  report.  User  documentation  and  the  final 
report  will  describe  the  actual  ISEIS  system  in  more  detail.  The  purpose  of  this  report  is  to 
describe  the  design  considerations  which  went  into  the  ISEIS  development  and  to  present  the 
theoretical  and  seismological  basis  for  the  processing  procedures  which  are  part  of  ISEIS. 

Section  1.0  gives  the  design  considerations  and  an  overall  top-level  design  of  the  system. 
Section  2.0,  taken  from  the  paper  by  Baumgardt  and  Young  (1990),  presents  the  results  of  a 
discrimination  study  of  western  Norway  mine  explosions  and  earthquakes.  The  results  of  this 
study  were  used  in  the  development  of  the  initial  event-identification  rules  of  ISEIS.  Section  3.0 
describes  the  Multiple  Event  Recognition  System,  or  MERSY,  which  we  have  developed 
specifically  for  the  identification  of  ripple-fired  explosions  using  the  observation  of  time- 
independent  spectral  modulations  or  spectral  "scalloping.”  Section  4.0  discusses  a  speech 
processing  algorithm,  called  dynamic  time  warping  (DTW),  which  we  have  adapted  for  seismic 
event  characterization  by  pattern  matching  of  regional  event  codas.  Finally,  Section  5.0  presents 
our  overall  conclusions  related  to  our  experience  in  the  design  of  ISEIS,  and  its  application  to  the 
characterization  and  identification  of  regional  seismic  events. 


SECTION  1.0 


DESIGN  CONSIDERATION  FOR  THE  INTELLIGENT  EVENT 
IDENTIFICATION  SYSTEM 


1.1  OBJECTIVES 


The  Intelligent  Seismic  Event  Identification  System  (1SE1S)  is  a  prototype  system  for  the 
systematic  identification  of  seismic  events  using  both  regional  and  teleseismic  seismic  data. 
Regional  data  from  regional  arrays  and  in-country  regional  networks  will  be  required  for  the 
identification  of  small  seismic  events.  However,  teleseismic  data  may  have  to  be  processed  for  the 
identification  of  larger  events  in  countries  for  which  in-country  seismic  network  data  is 
unavailable.  ISEIS  will  incorporate  the  latest  seismological  knowledge  of  seismic  discriminants, 
both  regional  and  teleseismic,  to  provide  as  complete  as  possible,  characterization  of  seismic  events 
for  identification.  In  addition  to  serving  as  the  primary  subsystem  for  event  identification  in  the 
Nuclear  Monitoring  Research  and  Development  (NMRD)  system,  ISEIS  is  also  being  designed  as 
an  analysis  tool  to  support  research  efforts  in  seismic  discrimination.  A  major  goal  is  to  develop  a 
system  which  is  flexible  enough  to  easily  incorporate  new  seismic  knowledge  and  event 
identification  techniques,  as  a  result  of  event  processing  and  event  identification  research. 

1.2  DESIGN  CONSIDERATION 


ISEIS  is  being  designed  as  a  post-processor  to  analyze  seismic  events  after  they  have  been 
detected  and  located,  using  some  combination  of  data  from  regional  arrays  and  seismic  networks  of 
single  stations.  We  assume  that  the  basic  event  parameters  (phase  identifications,  hypocenter. 
magnitudes)  have  been  determined  by  the  front-end  processors  of  the  NMRD  system  and  are 
available  in  the  ORACLE  relational  database.  The  task  of  ISEIS  is  to  take  as  input  the  event 
parameters,  to  make  additional  waveform  or  spectral  measurements  as  required  for  a  particular 


event-identification  method,  and  to  apply  knowledge  to  characterize  and  identify  the  event. 


The  basis  for  our  design  of  ISEIS  has  come  from  the  latest  results  in  seismic  discrimination 
research  as  well  as  knowledge  derived  from  other  system  development  projects.  Most  of  the  focus 
of  current  seismic  event  identification  has  been  in  regional  event  identification,  since  small-event 
identification  will  require  data  from  regional  stations  near  the  source.  Seismic  discrimination  using 
the  NORESS  regional  array  has  been  discussed  by  Pulli  and  Dysart  (1987),  Baumgardt  and  Ziegler 
(1989),  and  Baumgardt  and  Young  (1990).  Bennett  et  al  (1989)  has  presented  the  most  recent 
results  on  regional  event  identification  methods  using  regional  data  recorded  at  single  stations  for 
events  inside  the  Soviet  Union.  Regional  seismic  event  identification  is  uncertain  and  regional 
seismic  discriminants  are  highly  experimental,  which  means  that  event  identifications  using 
regional  data  must  be  qualified.  Baumgardt  (1990)  has  suggested  the  use  of  case-based  event 
characterization,  where  events  are  characterized  on  the  basis  of  how  similar  they  are  to  previous 
events  in  a  region.  This  differs  from  model-based  event  identification,  where  event 
characterization  in  terms  of  source  type  (i.e.,  explosion  or  earthquake)  is  made  based  on  a  source- 
type  model,  derived  either  from  first  principles  or  from  training  data.  Because  of  the  uncertainties 
of  regional  seismic-event  identification,  both  case-based  and  model-based  event  identification  will 
be  utilized. 

Teleseismic  discriminants  must  be  implemented  to  identify  events  which  could  occur  in 
third  countries  (i.e.,  outside  the  Soviet  Union).  Although  many  research  studies  have  tried  to 
apply  waveform  discriminants  (e.g.,  spectral  ratios,  complexity)  to  classify  seismic  events  as 
explosion  or  earthquake,  we  have  chosen  to  concentrate  primarily  on  physically^  based 
discriminants  which  are  reasonably  well  understood,  principally  location,  depth  and  the  Ms-mb . 

Finally,  our  design  is  based  on  results  of  ENSCO  efforts  as  part  of  the  Intelligent  Event 
Identification  (IAS)  project  (Bache,  1990).  This  includes  the  first  implementation  of  a  case-based 
reasoning  system,  called  seismic  script-matching  (Baumgardt,  1987;  Kandt  et  al,  1987),  and  the 
application  of  spectral  methods  for  identifying  ripple-firing  in  economic  explosions  to  distinguish 


them  from  earthquakes  (Baumgardf  and  Ziegler,  1989).  Also,  many  of  the  software  system 
development  concepts  of  the  IAS  development  project  have  been  applied  to  the  development  of 
ISEIS,  including  the  exploitation  of  distributed  systems  using  message  passing,  XI 1  graphics, 
database  techniques  using  SQL,  artificial  intelligence  methods,  and  rapid  prototyping. 

1.3  OVERALL  DESIGN  CONCEPT 

Figure  1  shows  the  overall  design  concept  currently  being  developed  in  ISEIS.  ISEIS  will 
provide  a  set  of  seismic-event  identification  techniques  or  discriminants  which  will  each  be  applied 
to  the  seismic  data  input  to  the  system.  As  shown  on  the  left  of  Figure  1,  a  standard  set  of 
"primary"  or  "active"  discriminants  can  be  applied  to  seismic  events  input  from  a  database  or 
selected  from  a  base  map  by  the  user.  Active  discriminants  are  those  which  are  judged  to  be 
reliable  and  relatively  well  understood.  Other  "secondary"  or  "inactive"  discriminants  will  also  be 
available  which  may  be  more  experimental  or  possibly  redundant  with  other  "active"  discriminants. 
"Inactive"  discriminants  can  be  "activated"  and  applied  to  an  event  by  the  user  in  an  interactive 
mode  for  research  purposes  or  to  provide  a  more  complete  characterization  of  the  events.  The 
discriminants  we  now  believe  to  be  the  most  robust  are  those  shown  under  Primary  Discriminants 
in  Figure  1.  (However,  research  by  Bennett  et  al  (1989)  and  Baumgardt  and  Young  (1990) 
indicate  that  the  L#  spectral  ratio  may  not  be  a  particularly  effective  discriminant  in  Scandinavia  and 
the  Soviet  Union,  althoughflt  has  proven  effective  in  other  regions.) 

In  addition  to  the  single  discriminants,  the  ISEIS  user  will  also  have  access  to  a  number  of 
signal-analysis  and  display  functions,  as  shown  on  the  right  of  Figure  1.  Some  of  these  functions 
w  ill  be  used  to  compute  features  needed  by  the  discriminants,  but  can  also  be  applied  by  the  user 
for  special  event  analysis  and  discrimination  research.  Results  of  these  processes  can  also  be  made 
available  to  the  event-identification  modules. 

The  center  column  shows  how  an  event  to  be  identified  will  flow  through  the  system.  The 
Event  Assessment  receives  as  inputs  event-parameter  and  waveform  data  information  from  the 
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FIGURE  I:  Intelligent  Seismic  Event  Identification 


Oracle  database.  The  system  will  access  this  information  from  Version  3.0  of  the  seismic  database 
of  the  Center  tor  Seismic  Studies  (CSS),  described  by  Anderson  et  al  ( 1990).  Kvent  parameters 
will  include  information  about  the  location  and  magnitudes,  stored  for  the  origin  and  origerr 
relations,  phase  identification  information  from  the  assoc,  and  arrival  relations,  and  information 
about  the  availability  of  waveform  data,  contained  in  the  w flag,  affiliation,  and  wfdisc  relations.  It 
i\  assumed  (hat  this  information  will  be  placed  in  the  database  by  the  front-end  detection, 
association,  and  location  processors  of  the  NMRD  system.  Event  Assessment  will  then  determine 
the  data  "status"  of  the  seismic  data  vis-a-vis  the  data  requirements  of  the  different  discriminants. 
The  status  will  be  determined  by  rules  in  an  expen  system  which  will  be  associated  with  each 
discriminant.  Status  rules  might  reason  on  such  things  as  the  minimum  number  of  stations 
required  for  a  discriminant,  whether  or  not  a  certain  type  of  phase  has  been  identified,  and  whether 
or  not  the  signal  to-noise  ratio  exceeds  a  certain  required  minimum,  etc.  The  status  assessment 
will  be  either  complete,  meaning  that  the  minimum  requirements  to  run  the  discriminant  are 
available,  incomplete  meaning  that  additional  data  is  required  to  run  the  discriminant,  and  no-data 
which  means  the  event  has  no  waveform  data  or  the  required  phase  identifications  have  not  been 
made.  An  explanation  of  the  status  will  then  be  presented  to  the  user  and  what  should  be  done  to 
upgrade  the  status,  if  necessary. 

In  Single  Discriminant  Analysis,  each  discriminant  is  applied  one  by  one  and  the  event  is 
characterized  using  either  case  based  reasoning  (CBR)  or  model-based  reasoning  (MBR).  or  both. 
In  MBR,  the  event  is  determined  to  be  either  earthquake  or  explosion  with  some  confidence  or 
classed  as  unknown.  In  CBR,  the  discriminant  feature  is  compared  to  those  for  reference  events 
and  the  event  is  classed  as  either  "typical"  or  "atypical"  of  the  reference  events,  with  an  associated 
confidence.  MBR  w  ill  be  accomplished  using  rules.  CBR  will  consist  of  matching  of  the  features 
of  the  new  event  to  templates  of  reference-events  and  the  confidence  of  match  will  be  computed, 
using  the  same  method  as  was  done  in  the  IAS  script  matcher,  in  Multidiscriminant  Analysis,  each 
discriminant  will  "vote"  on  the  identity  of  the  event,  using  both  MBR  and  CBR,  and  an  overall 
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assessment  of  the  events  characterization  will  be  made.  This  overall  characterization,  along  with 
the  specific  discrimination  results,  will  then  be  placed  in  the  database  as  feature  and  event- 
characterization  relations. 

An  expert  system  approach  will  be  used  to  reason  on  data  availability  for  status  assessment 
and  for  MBR  and  CBR  of  the  discriminant  processing.  The  expert  system  will  be  coded  in  the  C 
Language  Integrated  Production  System  (CLIPS),  an  expert  system  shell  developed  by  NASA  at 
the  Lyndon  B.  Johnson  Space  Center.  The  version  we  will  use  initially  will  be  Version  4.3, 
described  in  NASA  ( 1989).  An  evaluation  of  CLIPS  and  other  expen-system  shells  has  been  done 
by  Mettrey  (1991).  CLIPS  has  been  developed  in  the  C  language  and  is  designed  to  be  easily 
integrated  with  other  C  and  Fortran  languages.  In  ISEIS,  information  about  data  status  and  the 
results  of  the  discriminant  signal  processes  will  be  extracted  from  the  Oracle  database  using  C- 
language  SQL  queries.  This  data  will  then  be  passed  to  CLIPS  in  the  form  of  logical  "assertions" 
which  are  processed  by  the  CLIPS  rules  through  C  language  calls  to  routines  in  the  CLIPS 
libraries.  When  all  the  input  assertions  have  been  processed,  the  rules  will  make  conclusions 
which  will  be  output  in  the  form  of  textual  information  in  disk  files  that  will  later  be  displayed  to 
the  analyst.  In  addition,  the  CLIPS  rules  themselves  w  ill  call  other  C  language  routines  that  will 
write  the  results  to  special  relations  in  the  database.  The  text  files  and  the  database  relations  will 
then  be  accessed  by  the  ISL1S  explanation  functions. 

1.4  ACTOMATED  AND  INTERACTIVE  PROCESSING 

ISLIS  is  being  designed  to  be  run  in  either  an  automated  or  interactive  mode.  The 
functional  How  for  the  automated  processing  is  shown  in  Figure  2.  In  the  automated  mode,  a 
standard  set  of  acnve  discriminants  will  be  selected  which  will  be  run  on  all  the  events  placed  in  the 
database  by  the  front-end  NMRD  processors.  The  Database  Monitor  will  constantly  query  the 
database  for  new  events  which  have  been  added  and  must  be  identified,  subject  to  some  location 
and  size  constraints  (e  g.,  those  events  in  the  area  of  interest  with  mf,  or  Mi  greater  then  2.5).  The 
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FKJL’RE  2:  ISF.FS  Automated  Processing  Concept 


selected  events  are  then  passed  into  Technique  Data  Analysis  which  applies  the  active  set  of 
regional  and/or  teleseismic  discriminants.  These  discriminants  arc  applied  once  to  give  an  initial, 
quick-look  assessment  of  the  identity  of  the  event.  The  results  are  then  placed  into  the  database. 


Figure  3  shows  the  functional  flow  for  the  interactive  processing  in  ISEJS.  In  this  mode, 
events  formed  by  the  front-end  NMRD  processes  are  examined  by  the  analyst.  Events  of  interest 
will  also  be  plotted  in  the  ISTIS  Geographic  Display,  or  map.  along  with  other  geophysical  and 
seismological  information.  These  events  may  have  already  been  run  through  ISE1S  in  the 
automated  mode  and  may  have  initial  identifications.  The  analyst  can  then  select  events  for 
subsequent  processing,  which  are  again  passed  to  the  Tci  hnique  Data  Analysis  module.  As 
shown,  these  processes  can  be  applied  iteratively  or  repeatedly  to  further  refine  the  event 
characterization. 

In  the  interactive  mode,  the  events  will  then  be  presented  to  the  analyst  in  the  form  of  a 
status-results  spreadsheet  display.  A  highly  schematic  illustration  of  this  display,  as  it  will  appear 
in  X  1 1  on  a  SI  N  graphics  workstation,  is  shown  in  Figure  4.  The  data  status  and  discrimination 
results  are  displayed  to  the  user  in  the  form  of  numbers  and  color  coded  buttons  in  a  scrollable 
spreadsheet,  showing  a  matrix  of  discriminant  methods  across  the  lop  for  events  along  the  left- 
hand  column.  Each  square  in  Figure  4  provides  a  summary  display  for  the  status  and  results  for  a 
particular  discriminant,  shown  across  the  top.  applied  to  a  specific  event,  shown  plotted  down  the 
left  side.  Within  each  square  are  three  smaller  buttons  and  a  set  of  numbers.  The  left  button  gives 
discriminant  status  information,  which  refers  to  whether  or  not  the  data  required  by  the 
discriminant  o  available  to  the  system  The  middle  button  gives  MBR  results,  and  the  right  button 
provides  (  HR  results  The  numlvts  under  the  MBR  and  CBR  buttons  refer  to  the  event 
characterization  confidences  determined  by  the  techniques  for  the  event  Each  of  the  three  buttons 
is  active,  meaning  it  can  be  clicked  with  a  mouse  to  bring  up  other  displays  with  more  detailed 
information 
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FKIL'RK  4:  Schematic  Illustration  of  the  ISF2IS  Status-Results  Spreadsheet  Displat 


Each  of  the  buttons  in  Figure  4  will  be  color  coded  to  make  the  results  easier  for  the  analyst 
to  interpret.  With  reference  to  the  status  buttons  on  the  left,  if  the  button  is  green,  the  status  is 
"ready,"  meaning  all  data  is  available  and  the  discriminant  can  be  applied.  If  yellow,  the  status  is 
"incomplete."  meaning  only  part  of  the  data  required  by  the  discriminant  is  available.  The  button 
can  be  moused  to  provide  more  information  about  the  missing  data.  Finally,  if  the  button  is  red, 
then  no  data  is  available  and  the  discriminant  cannot  be  applied  to  the  event.  For  the  center  MBR 
button,  green,  red,  or  yellow  indicates  that  the  event  is  identified  as  either  an  earthquake, 
explosion,  or  unknown,  respectively.  For  the  CBR  button,  green,  red,  and  yellow  indicates  that 
the  event  is  typical,  atypical,  or  unknown,  respectively.  These  determinations  will  be  based  on  the 
confidence  of  match  of  the  discriminant  feature  to  reference  events  in  the  same  region  as  the  event. 
As  in  the  case  of  the  status  buttons,  the  MBR  and  CBR  buttons  can  be  moused  to  bring  up  more 
detailed  displays  of  the  processing  results  for  each  discriminant. 

The  spreadsheet  will  have  functions,  selectable  from  pull-down  menus  (not  shown  in 
Figure  4),  which  will  allow  the  analyst  to  apply  any  combination  of  discriminants  to  any 
combination  of  events,  or  to  apply  all  active  discriminants  to  all  events.  It  will  be  possible  for  the 
analyst  to  add  or  delete  events  or  discriminants  from  the  display.  Moreover,  the  analyst  can  invoke 
from  this  display  the  various  signal-analysis  functions  either  for  research  purposes  or  for  providing 
data  needed  by  the  discriminants. 

1.5  PROCESSING  REQUIREMENTS 

Our  design  of  the  ISFIS  system  has  been  based  on  our  assessment  of  the  current  state  of 
knowledge  in  seismic  event  identification  and  the  current  and  future  needs  of  a  monitoring  system 
for  a  Comprehensive  Testban  Treaty  (CTBT).  ISFIS  is  being  developed  as  an  integrated  system  to 
characterize  seismic  events  using  both  MBR  and  CBR  and  utilizing  both  regional  and  tele:  rismie 


Because  seismic-event  identification,  particularly  at  regional  distances,  is  not  well 
understood  and  is  a  field  of  active  research,  we  believe  it  is  important  that  a  case-based  reasoning 
capability  be  provided.  With  this  facility,  we  can  at  least  determine  if  the  event  looks  like  previous 
events,  even  if  a  discriminant  cannot  be  relied  on  to  determine  if  the  event  is  an  earthquake  or 
explosion. 

Although  regional  discrimination  will  be  of  greatest  importance  for  CTBT  monitoring. 
ISF.IS  may  still  be  required  to  process  events  from  third  countries  which  are  not  part  of  the  CTBT 
and  for  which  regional  seismic  data  may  be  unavailable.  Thus,  there  will  still  be  the  need  to  utilize 
teleseismic  data  for  these  regions.  Therefore,  ISKIS  will  have  the  capability  of  identifying  events 
using  teleseismic  data  for  larger  events.  Moreover,  for  larger  events,  combined  regional  and 
teleseismic  data  may  provide  more  definitive  event  characterizations  than  just  relying  only  on 
regional  or  teleseismic  data.  However,  in  the  initial  prototype,  the  primary  emphasis  will  be  placed 
on  the  identification  of  small  events  at  regional  distances. 

Finally,  the  system  we  are  developing  will  provide  a  useful  research  tool  forevent  analysis 
and  event  identification.  Our  modular  approach  to  the  development  and  the  spreadsheet 
implementation  ot  the  discriminants  will  allow  for  easy  modification  of  the  discrimination  methods 
as  new  knowledge  and  discriminant  techniques  are  developed  in  the  future. 

In  the  remainder  of  this  report,  we  discuss  in  detail  the  processing  requirements  for  ISFIS 
Section  20.  taken  from  Baumgardt  and  Young  (lWOi.  describes  the  seismological  basis  for 
regional  event  identification  m  western  Norway  using  regional  phase  amplitudes  and  spectral 
ratios.  The  initial  prototype  of  IS!  IS  will  be  designed  to  identify  events  in  western  Norway. 
Many  ot  the  ISFIS  display  s  will  tie  similai  to  those  show n  in  this  smily  Section  VO  describes  lhe 
sy  stem  which  will  he  included  to  identity  ripple  tired  explosions  Section  4  0  describes  a  new 
scheme  for  waveform  shape  matching  which  will  be  used  as  part  of  the  case  based  reasoning 
process  Section  S  O  presents  some  final  conclusions  relevant  to  the  design  of  ISFIS 


SECTION  2.0 


REGIONAL  WAVEFORM  DISCRIMINATION 
2.1  INTRODUCTION 

Effective  seismic  event  identification  of  nuclear  explosions  and  earthquakes  has  long  been  a 
ma  jor  goal  in  testban  treaty  monitoring,  as  evidenced  by  the  many  studies  reviewed  by  Pomeroy  et 
al  ( 1982)  which  were  done  up  to  the  early  19X0s.  Many  subsequent  studies  have  been  undertaken, 
motivated  to  a  large  extent  by  the  increased  availability  of  regional  seismic  data,  including  regional 
array  data.  In  spite  of  all  the  research  in  this  area,  no  consistently  reliable  and  universally 
applicable  regional-waveform  discriminant  has  yet  been  discovered. 

Part  of  the  problem  has  been  a  lack  of  a  unifying  theory  of  regional  seismic  event 
identification  which  completely  explains  how  seismic  source  phenomenology  affects  recorded 
seismic  signals,  l  or  example,  simple  source  physics  of  earthquakes  and  explosions  would 
suggest  that  earthquakes,  being  dislocation  sources,  should  generate  more  shear-wave  energy 
relative  to  compressional-wave  energy  than  explosions,  and  many  of  the  studies  reviewed  by 
Pomeroy  et  al  (19X2)  that  tested  S-to-Pratio  discriminants  showed  some  separation  between 
earthquakes  and  explosions.  However,  studies  of  the  P/Eg  ratio  discriminant  by  N'uttli  (1981), 
Bennett  and  Murphy  (19X6).  and  Taylor  et  al  (19X9)  have  found  significant  overlap  in  the 
earthquake  and  explosion  populations  for  events  in  Eurasia  and  in  western  U  S.  However,  more 
recent  studies  by  Pulli  and  Dysart  ( 19X7),  for  events  in  Scandinavia,  and  Bennett  et  al  ( 1989),  for 
events  in  the  Eurasian  shield,  have  shown  P/Eg  type  ratios  to  be  effective  for  separating 
earthquakes  and  explosions  Evidently,  differences  in  propagation  paths  may  have  caused  the  poor 
performance  of  this  discriminant  in  some  of  the  studies,  but  the  exact  nature  of  the  trade-off 
between  source  and  propagation  path  effects  on  the  P/Eg  ratio  discriminant  is  not  well  understood. 

Mixed  results  have  also  been  reported  for  spectral  discriminants.  Murphy  and  Bennett 
( 19X2)  and  Bennett  and  Murphy  ( 19X6)  found  that  NTS  explosions  and  western  U.S.  earthquakes 


13 


separated  on  Eg  spectral  ratios,  with  earthquake  Li>  waves  having  more  high  frequency  content 
than  explosions.  Similar  results  were  obtained  by  Taylor  et  al  ( 19XX)  for  Eg  and  other  phases  for 
events  in  tne  same  region.  However,  Pulli  and  Dysart  ( 19X7)  and  Bennett  et  al  (19X9)  reported 
that  the  discriminant  is  less  effective  in  separating  earthquakes  and  explosions  in  Scandinavia  and 
in  the  eastern  European  shield,  respectively.  Again,  as  pointed  out  by  Bennett  et  al  (19X9). 
propagation  path  differences  for  the  explosion  and  earthquake  populations  may  have  an  effect  on 
the  performance  of  this  discriminant. 

Another  problem  associated  with  small-event  identification  is  distinguishing  economic 
blasting  from  nuclear  explosions  and  earthquakes.  Baumgardt  and  Ziegler  (19XX)  and  Hedlin  et  al 
(  19X9,  1990)  have  shown  that  economic  blasting  can  be  identified  by  observing  persistent  spectral 
modulations  produced  by  ripple-fire.  However,  this  discriminant  could  be  spoofed  since  nuclear 
explosions  can  also  be  ripple-fired,  although  it  might  be  difficult.  For  example,  Baumgardt  and 
Ziegler  ( 19XX)  found  that  NORSAR  recordings  of  presumed  peaceful  nuclear  explosions  in  Eurasia 
had  the  same  persistent  modulations  observed  in  ripple-fired  mine  blasts,  indicating  that  the  nuclear 
explosions  were  ripple-fired.  Moreover,  Bennett  et  al  ( 19X9)  argue  that  spectral  modulations  are 
unobservable  for  known  mine  blasts  in  the  L'.S..  perhaps  because  the  ripple-fire  delay  times  were 
too  short  to  be  observed  in  the  limited  bandwidth  of  the  data.  The  observation  of  spectral 
mixlulations  may  identify  many,  perhaps  most,  economic  explosions,  depending  on  number  and 
delay  times  of  the  ripple-fired  explosions  and  the  bandwidth  of  the  recording  instrumentation. 
However,  there  are  conceivable  scenarios  where  the  discriminant  may  fail. 

Taken  together,  these  studies  have  shown  that  regional  discriminants  cannot  be  applied  in 
the  same  wav  everywhere,  that  thev  are  highly  dependent  on  the  nature  of  the  regional  phase 
propagation  path  effects,  and  that  there  is  a  strong  regional  variability  in  the  effectiveness  of 
discriminants.  Because  of  this,  a  case  based  approach  may  have  to  be  applied,  where  the 
waveform  characteristics  of  an  event  are  compared  with  those  of  previously  observed,  known 
events  and  identified  on  the  basis  of  the  comparison,  assuming  that  the  propagation  paths  arc- 
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common.  Ideally,  the  unknown  event  and  the  known  reference  events  must  be  in  the  same  region 
so  that  propagation-path  differences  do  not  bias  the  discrimination  results.  However,  it  has  been 
difficult  in  previous  discrimination  research  to  find  populations  of  earthquakes  and  explosions 
which  are  collocated  and  there  is  no  guarantee  that  in  practice  collocated  reference  events  of  known 
identity  will  be  available  to  compare  with  unknown  events.  Thus,  in  an  operational  event 
identification  scenario,  the  seismic  analyst  will  be  forced  to  work  with  what  is  available  and  trv  to 
account  for  an\  propagation  effects  which  might  bias  event  identification  using  previously 
observed  cases. 

In  this  paper,  we  w  ill  explore  some  of  the  problems  associated  with  using  the  case-based 
approach  to  event  identification.  Baumgardt  ( 1990)  has  discussed  the  application  of  an  artificial 
intelligence  technique,  known  as  "case  based  reasoning."  for  systematically  identifying  seismic 
events  on  the  basis  of  comparison  with  cases.  In  this  paper,  we  will  focus  primarily  on  the 
seismological  and  signal  processing  aspects  of  the  problem.  Our  emphasis  will  be  on  the  analysis 
of  regional-array  data  from  the  NORKSS  array  and  the  discrimination  of  earthquakes  and  mine 
explosions  located  in  western  Norway.  We  will  discuss  the  question  of  whether  propagation-path 
effects  bias  waveform  discriminants  and  what  considerations  need  to  be  made  when  attempting  to 
identify  actual  case  events. 

2.2  DATA  AND  PROCESSING  METHODS 
\ORESS  Database 

1  or  this  study.  NORKSS  data  for  a  group  of  earthquakes  and  economic  explosions  located 
in  western  Norway  were  analyzed,  figure  1  shows  the  locations  of  the  presumed  earthquakes, 
labeled  as  Q,  and  the  blasting  sites.  Blasjo  (Bl.A).  a  dam  excavation  site,  and  the  Titania  (TITA) 
mine.  The  source  parameters  of  the  earthquakes  and  blasts  are  shown  in  Tables  1  and  2. 
respectively.  All  the  events  in  Table  1  were  reported  in  the  Bergen  regional  bulletin,  except  Q12 
which  was  presumably  below  the  detection  threshold  of  the  Bergen  network.  The  location  of  Q12 


in  Table  1  was  determined  by  NORESS  The  basis  for  the  presumption  of  earthquakes  for  these 
events  is  simply  that  they  were  not  reported  as  explosions  in  the  Bergen  bulletin  or  otherwise 
known  to  be  blasts  at  active  mines.  All  the  events  in  Table  2  were  reported  in  the  Bergen  bulletins 
and  confirmed  to  be  blasts  at  the  Titania  mine  and  the  Blasjo  dam  excavation  site. 

In  order  to  examine  the  effect  of  event  location  on  seismic  waveform  features,  the  events  in 
Table  I  and  in  Figure  5  have  been  divided  into  five  regions.  The  Region  1  events  occurred  off¬ 
shore  near  the  MaloyT.'lstein  region  Event  Q1  has  the  highest  local  magnitude  reported  in  the 
Bergen  bulletin  (4.2)  and  was  also  reported  in  the  Preliminary  Determination  of  Epicenters  bulletin 
as  having  a  body  wave  magnitude  of  5.0.  This  event  was  one  ot  the  largest  events  to  occur  in 
Norway  in  the  past  40  years  and  was  felt  over  most  of  southern  and  central  Norway  (Hansen  et  al, 
1980).  The  events  which  occurred  shortly  afterwards  (Q2-Qb)  were  aftershocks  of  Ql,  whose 
spectral  scaling  properties  were  studied  by  Chael  and  Cromer  ( 1988). 

The  Region  2  events  (Q9.  Ql(),  and  Ql  1.)  ail  occurred  on-shore  in  the  Maloy  region.  The 
question  marks  indicate  that  we  suspect  they  may  actually  be  explosions  rather  than  earthquakes, 
for  reasons  that  will  be  discussed  later  Region  4  has  one  event  located  to  the  northeast  of  the 
events  in  Region  2.  T  his  small  event  was  detected  at  NORESS.  but  not  by  the  Bergen  network, 
and  the  source  parameters  in  Table  I  were  determined  by  NORESS.  We  also  believe  this  event  is  a 
blast  rather  than  an  earthquake. 

Region  4  contains  two  events  (Q14  and  Q15)  which  occurred  in  the  vicinity  of  Bergen  and 
a  third  event  (Ql  4>  which  occurred  further  sou’h  in  the  Stavanger  region.  The  Q14  event,  with  the 
largest  magnitude  ot  2.9.  may  have  been  felt  locally  (Erode  Ringdal  and  Svein  Mykkeltveit. 
personal  communication).  These  are  the  closest  events  to  the  TIT  A  and  BEA  sites  and  were  not 
reported  as  blasts,  although,  as  shall  be  discussed  later,  we  suspect  that  the  Stavanger  event.  Q14. 
mav  actual iv  be  a  blast 


lb 


NORESS 


FKiliRK  5:  Map  showing  the  locations  of  the  presumed  earthquakes,  labelled  as 
Q,  and  the  locations  of  the  Rlasjo  and  Titania  blast  sites.  The  brackets  refer  to 
the  regionalization  of  the  events  discussed  in  the  text.  All  event  locations,  except 
for  the  event  in  Region  3  (QI2),  were  determined  by  the  regional  seismic  network 
of  the  I  niversity  of  Bergen.  The  location  of  QI2  comes  from  NORESS. 
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TABLE  1 


EPICENTERS  OF  PRESUMED  WESTERN  NORWAY  EARTHQUAKES 
AND  OTHER  UNKNOWN  EVENTS 


Event 

Date 
(m/d/y ) 

Origin  Time 
(UTC) 

Latitude 

(deg.N) 

Longitude 
(deg.  E) 

Ml 

Dist. 

(km) 

Region  l 

Q1 

2/05/86 

17:53:16.2 

62.74 

4.63 

4.2 

429 

Q2 

2/05/86 

18:50:03.4 

62.27 

4.69 

2.8 

403 

Q3 

2/05/86 

20:23:29.8 

62.41 

6.06 

2.7 

303 

Q4 

2/05/86 

20:31:37.0 

62.79 

4.59 

2.2 

433 

Q5 

2/05/86 

23:35:41.1 

62.74 

4.50 

2.6 

434 

Q6 

2/06/86 

06:19:52.4 

62.90 

4.86 

2.3 

427 

Q7 

02/13/86 

13:39:00.3 

62.40 

5.28 

2.5 

381 

Q8 

02/13/86 

19:03:48.2 

62.61 

5.07 

2.6 

401 

Region  2 

Q9(?) 

2/05/86 

15:57:02.4 

62.05 

5.37 

2.0 

361 

Q10(?) 

02/16/86 

18 

19:41.3 

61.69 

4.90 

2.0 

373 

Q1K?) 

02/14/86 

16:51:05.1 

61.68 

4.97 

1.8 

369 

Region  .1 

Q12(?) 

2/05/86 

15 

22:44.0 

62.5 

6.82 

1.6 

325 

Region  4 

Q13  r>) 

12/07/85 

14:39:09.0 

58.90 

5.98 

2.0 

373 

Q14 

1 1/27/85 

04:53:32.8 

59.73 

5.71 

2.9 

342 

Q15 

02/15/86 

18:31:46.4 

59.86 

5.73 

2.1 

336 

Region  5 

UND1 

1 1/20/85 

22 

10:44.2 

57.61 

5.67 

2.3 

483 

UND2 

1  1/20/85 

22 

24:38.1 

57.66 

5.72 

2 _2 

478 

UND3 

1 1/20/85 

22 

57:10.8 

57.63 

6.27 

2.2 

459 

UND4 

1 1/20/85 

23 

10:47.5 

57.66 

5.35 

2.3 

493 

UND5 

1 1/20/85 

23 

17:28.9 

57.69 

5.45 

2.3 

486 

l\D6 

1 1/20/85 

23 

23:10.0 

57.64 

5.62 

2.2 

483 
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TABLE  2 


EPICENTERS  FOR  CONFIRMED  ECONOMIC  EXPLOSIONS 


Event  Date  Origin  Time  Size  Ml 


(Blasjo  Explosions  59. J 1°  N  6.95°  E  Dist.=  JO  I  km) 


f:x  i 

08/05/83 

17:42:58.7 

62.9 

2.6 

f;x2 

08/06/85 

17:50:07.9 

30.8 

2.4 

KX3 

10/17/85 

10:00:00.4 

32.7 

2.4 

(Titania  Mine 

Explosions  5X342°  N  6.425°  E  Dist. 

=  J94 km) 

HX4 

1 1/08/85 

14:18:54.6 

132.5 

2.4 

F.X5 

02/14/86 

14:13:24.9 

95.7 

2.7 

HX6 

02/14/86 

17:54:10.6 

16.2 

2.3 

KX7 

01/07/86 

14:14:28.9 

43.5 

2.2 

HX8 

01/17/86 

14:11:01.5 

43.9 

2.7 

EX9 

01/07/88 

14:24:43.5 

77.4 

2.2 

EX  10 

02/10/88 

14:17:46.5 

103.2 

2.5 

EX  11 

03/17/88 

14:13:10.3 

95.1 

2.4 

EX  12 

03/28/88 

13:17:27.0 

74.2 

2.4 

Finally,  we  have  included  a  group  of  events  in  Region  5,  near  the  TITA  blast  site,  which 
we  have  indexed  in  Table  1  as  UND1  through  UND6.  Suteau-Henson  and  Bache  (19X8)  studied 
these  events  and  compared  Lg  spectral  ratios  of  these  events  with  those  of  events  at  TITA.  They 
suggested  that  these  events  were  earthquakes  because  they  occurred  off-shore.  However,  we  shall 
show  data  later  which  suggests  that  these  events  may  have  actually  been  underwater  explosions, 
which  is  why  we  refer  to  them  by  the  UND#  index. 

Incoherent  Beam  Analysis 

Incoherent  beams  on  bandpass  prefiltered  waveforms  were  used  to  measure  regional-phase 
amplitudes.  Incoherent  beamfonming  consists  of  computing  log-rms  amplitudes  in  adjacent,  one 
second  time  windows  on  each  channel  of  the  array,  starting  about  one  minute  before  the  Pn- wave 
onset  time  and  extending  through  the  seismogram  into  the  Ly  coda.  The  log-rms  amplitudes  for 
each  time  window  are  then  averaged  across  all  the  array  elements.  When  plotted  versus  time, 
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incoherent  beams  give  an  envelope  description  of  the  seismic  trace  viewed  in  log-amplitude,  or 
magnitude,  space. 

Figures  6  and  7  show  filtered  waveforms  from  the  NRAO  element  of  NORESS  and 
incoherent  beams,  computed  using  all  vertical,  short  period  NORESS  elements,  plotted  on  the 
same  time  axes,  for  an  earthquake  and  explosion.  The  horizontal  dashed  lines  on  the  incoherent 
beams  are  the  average  rms  noise  levels  over  a  one  minute  time  interval  ahead  of  the  Pn  onset.  The 
8-16  Hz  incoherent  beam  and  average  noise  level  have  been  shifted  up  for  visibility  relative  to  the 
2-4  Hz  beam.  The  presumed  onsets  of  the  regional  phases,  Pn,  Rg,  Sn,  and  Eg,  are  indicated  on 
both  the  waveform  and  incoherent  beam  plots.  Comparison  of  the  waveform  plots  with  the 
incoherent  beams  indicates  that  some  phases,  notably  Sn,  are  easier  to  see  on  the  incoherent  beam 
plots.  A  comparison  of  Figures  6  and  3  reveals  that  the  Pn  wave  is  not  visible  in  the  low 
frequency  (2  to  4  Hz)  filter  band  for  the  earthquake  in  Figure  6  whereas  it  is  apparent  in  both  the  2 
to  4  Hz  and  8  to  16  Hz  band  for  the  explosion  in  Figure  7.  This  observation  was  first  pointed  out 
by  Baumgardt  and  Ziegler  (1988)  and  has  important  implications  for  discriminating  these  events, 
which  will  be  discussed  in  the  next  section. 

Spectral  A  nalysis 

Our  spectral  analysis  method  is  the  same  array  stacking  procedure  described  by  that  of 
Baumgardt  and  Ziegler  (1988).  In  brief  the  Fourier  power  spectrum  is  computed  on  each  channel 
on  the  windowed  phases  and  on  the  noise  background  to  Pn.  The  power  spectra  on  each  channel 
are  then  corrected  for  noise  and  instrument  and  averaged  across  the  array.  The  window  lengths  for 
each  phase  varied  depending  on  the  duration  of  each  phase,  but  in  general  were  between  7  and  14 
seconds  for  Pn,  P%.  and  Sn,  and  25.6  seconds  for  Eg  and  noise. 
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Incoherent  Beams 


Pn  Pg  Sn  Lg 


Time  CsecD 


FIGURE  6:  Top:  Waveforms  recorded  at  the  NRAO  array  element  of  NORESS 
for  the  Q4  earthquakes  after  bandpass  filtering  in  the  2-4  and  8-16  Hz  bands. 
Bottom:  Incoherent  beams  using  one  second  averaging  windows,  computed  from 
all  25  NORESS  vertical  component  traces  after  prefiltering  the  traces  in  the  2-4 
and  8-16  Hz  bands. 
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2.3  DISCRIMINATION  ANALYSIS  RESULTS 


Regional  Phase  Spectra 

Figures  8,  9,  and  10  show  selected  array-averaged  spectra  for  the  regional  phases,  Pn ,  Pg, 
Sn,  and  Lg,  for  events  in  the  five  regions.  In  each  case,  the  spectra  for  the  Pg,  Sn,  and  Lg  phases 
were  shifted  upward  by  0.5  log  units  relative  to  the  Pn  and  noise  spectra  for  purposes  of  visibility. 
The  spectra  of  the  Pn  background  noise  are  plotted  as  dashed  lines. 

Figure  8a  shows  the  spectra  for  one  of  the  Region  1  events  compared  with  spectra  for  an 
explosion  in  Figure  8b.  All  the  spectra  in  Region  1  resembled  that  in  Figure  8a  in  that  the  regional 
phase  spectra  in  general  were  very  simple  and,  at  frequencies  above  2  Hz,  above  the  frequencies 
where  the  low-frequency  effects  of  the  instrument  removal  are  apparent,  the  spectra  fall  off  linearly 
with  frequency  into  the  noise.  In  contrast,  the  explosion  spectra  in  Figure  8b  exhibit  strong 
modulations  or  scalloping,  which  Baumgard t  and  Ziegler  (1988)  attributed  to  multiple  explosions 
or  "ripple-fire."  The  key  feature,  which  is  diagnostic  of  ripple-fire,  is  that  the  same  modulation  is 
apparent  in  all  spectra,  hence,  the  modulations  are  time  independent.  Hedlin  et  al  (1989)  have 
shown  that  the  time-independent  modulations  are  apparent  in  coda  waves  as  well  as  in  regional- 
phase  spectra.  Most  of  the  events  that  we  have  studied  at  TITA  and  BLA  showed  obvious  spectral 
modulations.  The  exception  was  the  January  17,  1985  event  (EX  10),  which  showed  no  evidence 
of  modulations.  Thus,  modulations  may  or  may  not  be  present  in  blasts,  and  when  present,  may 
vary  in  intensity  and  periodicity,  depending  on  how  the  ripple-fire  pattern  is  designed  and  the 
bandwidth  of  the  recording  seismometers. 

Figure  9a  shows  the  spectra  of  one  of  the  three  events  in  Region  2  (Q9).  This  event  show  s 
a  arong  modulation  pattern  apparent  in  all  spectra  which  could  only  have  been  produced  by  a 
multiple  event.  Event  Q 1 2  in  Region  3  had  a  very  similar  set  of  spectra,  whose  modulations 
closely  resemble  those  observed  for  known  blasts.  Figure  9b  shows  a  second  event  in  Region  2 


(b) 


FKil  RK  8:  (a)  Regional  phase  spectra  for  a  regional  event  (Q3)  in  Region  1, 

strongly  suspected  to  be  an  earthquake.  The  Pg,  Sn  and  Lg  spectra  have  been 
shifted  up  by  0.5  log  amplitude  units  relative  to  the  noise  and  Pn  spectrum,  (b) 
Regional  phase  spectra  for  a  confirmed  explosion  (EX9)  located  at  the  Titania 
mine,  which  exhibits  time-independent  spectral  modulations  indicative  of  ripple¬ 
firing. 
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FKil'RK  9:  (a)  Regional  phase  spectra  for  a  regional  event  (Q9)  in  Region  2. 

Originally  presumed  to  he  an  earthquake,  this  event  has  time-independent  spectral 
modulations  indicative  of  ripple-fire,  (b)  Regional  phase  spectra  for  a  regional 
event  (QIO)  in  Region  2,  which  has  an  indication  of  a  single  half  cycle  of  a 
modulation  pattern  for  a  ripple-fire  with  delay  time  less  than  0.05  seconds. 
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FUJI  RK  10:  la)  Regional  phase  spectra  for  the  southern-most  event  in  Region  4 
iQ  1 3 1 .  Originally  presumed  to  be  an  earthquake,  the  event  exhibits  time- 
independent  spectral  modulations  indicative  of  ripple-fire,  (b)  Regional  phase 
spectra  for  an  e\ent  in  Region  5  < LJND5),  which  appears  to  be  an  underwater 
explosion.  The  time-independent  modulations  may  be  due  to  bubble-pulse 
interference  and/or  underwater  reverberations. 
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(Q10)  which  has  a  broad  hump  from  about  3  H/  to  near  20  M/.  This  may  be  a  very  subtle 
modulation  produced  by  multiple  events  delayed  by  times  near  or  less  than  0.05  seconds,  which 
would  produce  null  near  or  beyond  the  \yi|uist  frequency  of  20  11/  F.vent  Qll  had  regional 
phase  spectra  similar  to  those  of  ()  I  I 

f  igure  10a  shows  the  spectra  of  the  Stavanger  event  Q 1  3.  These  spectra  have  a  strong  null 
at  between  10  and  1  1  H/.  which  appears  to  be  caused  by  a  spectral  modulation.  Figure  10b  shows 
one  of  the  spectra  of  the  off-shore  events  near  TITA  (l  \D5)  1'his  event  and  all  the  L’N'D  events 
had  very  strong  spectral  modulations.  Interestingly,  the  periodicity  of  the  modulations  for  all  the 
events  in  this  group  were  nearly  the  same. 

The  key  question  which  must  be  addressed  when  considering  whether  time-independent 
spectral  modulations  imply  ripple-fired  blasts  is  whether  or  not  earthquakes  can  produce  time- 
independent  modulations  as  well.  In  principle,  there  is  nothing  which  precludes  earthquakes 
consisting  of  multiple  ruptures,  which  is  commonly  observed  for  large  earthquakes.  Blandford 
( 1*475).  for  example  has  modeled  complex  earthquakes  as  consisting  of  the  superposition  of  many 
subearthquakes "  in  order  to  explain  a  number  of  teleseismic  discriminants.  It  is  unknown 
whether  such  a  model  might  also  apply  to  small  earthquakes  and  whether  such  complex 
earthquakes  can  produce  the  kinds  of  coherent  time  independent  modulations  that  have  been 
observed  for  ripple-fired  explosions  Thus,  the  presence  of  time  independent  spectral  modulations 
alone  does  not  prove  that  events  Q*>.  UI0,  Ql  l.  (<>12.  and  QI3  are  ripple  fired  blasts. 

It  is  more  likely  that  1  \l)l  through  1  M)0  are  m  tact  underwater  blasts.  First,  the  origin 
times  ot  the  events,  although  uncertainly  determined  by  the  Bergen  network,  are  very'  evenlv 
spaced  in  time  by  about  14  minutes  Second,  the  time  independent  modulations  in  all  the  spectra 
indicate  that  the  multiple  events  were  delayed  bv  nearly  the  same  time,  which  is  about  0.3  seconds. 
It  is  unlikely  that  eitliet  ot  these  two  phenomena  would  be  observed  in  earthquakes.  Baumgardt 
and  Ziegler  t  have  suggested  that  the  spectral  modulations  in  these  events,  if  they  are 


underwater  explosions,  may  have  been  caused  by  interference  of  bubble  pulses  in  the  water  or 
reverberations  in  the  acoustic  waveguide  produced  by  the  water  column.  The  consistency  of  the 
rime-independent  modulations  indicates  that  each  of  the  underwater  blasts  were  probably  detonated 
at  about  the  same  depth  in  the  water. 

Comparison  of  Incoherent  Beam  Amplitude  Ratios 

As  has  been  discussed  above,  the  relative  excitation  of  compressional  and  shear  wave 
energy,  represented  in  terms  of  /Mo -S  ratios,  has  been  considered  as  a  possible  discriminant 
between  explosions  and  earthquakes.  In  theory,  earthquakes  should  generate  more  shear  energy 
relative  to  compressional  energy  than  explosions.  However,  propagation  effects  must  also  be 
considered  when  comparing  these  ratios  for  populations  of  earthquakes  and  explosions.  In  this 
section,  we  examine  regional  phase  amplitude  ratios  for  NORESS  recordings  of  blasts  and 
earthquakes  on  a  region-by-region  basis  in  order  to  consider  possible  regional  variations  in  these 
features. 

Frequency-dependent  amplitude  ratios  between  the  phase  pairs,  Pn  and  Pg,  Pn  and  Sn,  Pn 
and  /.g.  Pn  and  Sn.  and  Pg  and  /.g.  were  determined  from  the  incoherent  beams,  discussed  above. 
First,  incoherent  beams  were  computed  for  the  vertical  component  traces  after  prefiltenng  the 
seismograms  using  a  set  of  six-pole.  Butterworth  recursive  filters.  The  prefilter  bandpasses  were 
2.0  to  4.0.  2.5  to  4.5,  3.0  to  5  0.  4.0  to  6.0.  5.0  to  7.0,  6.0  to  S  O.  S.O  to  10.0  and  8.0  to  16.0 
11/  The  beam  traces  were  then  plotted  and  the  times  of  the  Pn ,  Pg.  Sn.  and  /.g  peaks  were  noted, 
as  shown  in  Figures  6  and  7.  We  have  found  all  these  phases  to  be  most  distinct  on  the  incoherent 
beam  in  the  X  to  16  11/  filter  band.  Therefore,  we  picked  the  peak  amplitudes  in  this  filter  band  and 
used  the  same  rimes  to  measure  the  amplitudes  in  the  other  filter  bands. 

For  these  events,  we  have  two  choices  for  regional  P.  Pn  and  Pg.  and  regional  S.  Sn  and 
I  n.  tor  computing  amplitude  ratios  Baumgardt  ( 1  W0>  considered  the  frequency  dependence  of  all 
possible  ratios.  Pn  Pn.  Pn  Sn.  Pn  l. g.  Pg/.Vn.  and  Pg/Lg,  and  found  that  the  greatest  difference 


between  the  blast  and  earthquake  groups  was  apparent  in  the  Pn  Sn  and  Pn'L g  ratios  at  high 
frequency  (K  16  H/»  f  igures  I  1  and  12  show  the  Pn  Sn  and  Pnil.g  ratios,  respeetively.  on 
ineoherent  beams  tor  the  X- 161  1/  pret'ilter  plotted  separately  tor  each  of  the  five  regions  and  for  the 
two  blast  sites.  Note  that  the  Bi.A  blast  site  is  near  Region  4  and  the  T1TA  blast  site  is  close  to 
Region  5. 

Different  symbols  are  also  used  to  indicate  the  different  kinds  of  sources.  The  triangles 
indicate  the  events  we  are  reasonable  sure  are  earthquakes  and  did  not  exhibit  time-independent 
spectral  modulations.  The  diamonds  indicate  even's  which  were  originally  thought  to  be 
earthquakes,  but  whose  spectra  appear  to  have  time-independent  spectral  rrodu’  io'  which 
resemble  those  produced  by  ripple  tired  explosions.  The  circles  represent  all  the  confirmed  blasts, 
most  of  which  had  time- independent  spectral  modulations. 

Both  f  igures  1  1  and  12  show  that  the  events  in  Region  1  and  two  events  in  Region  4.  all  of 
which  we  strongly  suspect  tire  earthquakes,  have  very  low  values  of  Pn/Sn  and  Pn/Lg  ratios.  It  is 
interesting  to  note  that  the  variance  of  the  ratios  in  Region  1  is  very  low,  in  spite  of  the  fact  that  the 
Bergen  locations  of  these  events  are  distributed  out  over  a  30  to  50  km  area.  Also,  the  amplitude 
ratios  for  the  two  Region  4  presumed  earthquakes,  about  300  km  south  of  the  Region  1  events,  are 
the  same  as  the  Region  1  earthquakes.  Thus,  for  the  events  which  are  most  likely  earthquakes,  the 
amplitude  ratios  are  consistent  and  do  not  seem  to  depend  on  differences  in  propagation  path  from 
the  two  regions  to  NORTSS. 

The  confirmed  mine  blasts  have  much  higher  ratios  than  the  earthquakes.  Also,  the 
variance  in  these  estimates  for  the  same  mine  are  much  larger  than  those  of  the  earthquake  group, 
even  though  the  blasts  at  each  site  supposedly  occurred  at  the  same  location.  The  BLA  blasts  also 
appear  to  have  greater  excitation  of  Pn  relative  to  Sn  and  l.g  than  the  TITA  events.  However,  the 
confirmed  blasts  appear  to  K-  clearly  separated  (rom  the  earthquake  group,  with  the  blasts  having 
greater  compression.!)  wave  energy  relative  to  shear-wave  energy  than  the  earthquakes. 
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H(il  Rh  II:  Plot  of  the  values  of  the  Pn/Sn  ratio  in  the  8-16  H/,  frequency  band 
for  each  of  the  five  regions  and  for  the  two  confirmed  economic  blast  regions, 
BLA  and  TIT  A.  The  values  of  Pn  and  Sn  amplitudes  were  measured  from 
incoherent  beams  of  prefiltered  waveforms. 
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H(.'  K!'  12:  Same  us  Figure  II  for  the  Pn  l.g  amplitude  ratios. 


It  is  also  apparent  in  Figures  1 1  and  12  that  the  events  in  Regions  2,  3,  and  4,  which  we 
suspected  as  being  explosions  based  on  their  having  time-independent  modulations  in  the  spectra, 
have  higher  PniSn  and  PnlLg  r  os  than  the  earthquakes,  although  they  are  not  quite  as  high  as  the 
average  of  all  the  blast  ratios.  However,  they  do  overlap  the  bottom  pan  of  the  TIT  A  population, 
as  might  be  expected  since  the  distances  of  the  events  in  these  regions  are  close  to  those  of  the 
TITA  blasts. 

The  suspected  underwater  explosions  in  Region  5  have  higher  ratios  than  the  earthquakes, 
and  they  are  about  the  same  as  those  of  the  BLA  events.  Also,  the  variance  in  the  amplitude  ratios 
for  the  Region  5  group  is  much  smaller  than  those  of  either  the  TITA  or  BLA  groups.  The  higher 
excitation  of  Pn  relative  to  Sn  and  Lg  compared  with  the  other  events  is  consistent  with  the  UND 
events  being  blasts  in  the  water,  since  these  events  should  have  generated  no  intrinsic  shear-wave 
energy.  Thus,  the  Sn  and  Lg  energy  comes  entirely  from  mode  conversions  at  the  water-bottom 
interface  w  ith  the  solid  earth  and  by  scattering  in  the  earth. 

In  order  to  determine  how  the  amplitude  ratios  depend  on  the  absolute  levels  of  Pn ,  Sn,  and 
Lg.  Figure  13  shows  scatter  plots  of  the  logs  of  the  ratios  versus  the  logs  of  the  absolute  rms 
amplitudes.  For  both  the  PniSn  and  PnlLg  ratios,  we  find  that  the  variations  in  the  ratios  seem  to 
be  slightly  more  correlated  with  the  Pn  amplitudes  than  with  the  Sn  or  Lg  amplitudes,  particularly 
within  the  blast  and  modulated-spectra  group.  Notice  also  that  all  the  events,  save  one,  which  is 
the  large  foreshock  (Q1 ),  have  nearly  the  same  Sn  and  Lg  amplitudes  and  that  the  earthquake 
population  cleanly  separates  from  the  explosion  and  modulated  spectra  group  on  the  basis  of  the 
Pn  Sn  and  Pml.g  ratios.  This  was  expected  since  the  events  were  selected  for  a  limited  local 
magnitude  range,  based  on  the  Bergen  coda  duration  magnitude  measure.  Thus,  we  conclude  that 
most  of  the  variation  in  the  ratios  comes  from  variations  in  the  Pn  excitation,  not  the  Lg  excitation. 

Our  results  are  consistent  with  those  of  other  studies  that  the  Pn  wave  is  stronger,  relative 
to  the  other  phases,  for  explosions  than  for  earthquakes  and  that  the  discriminatory  capability  of  Pn 


32 


Log  Pn/Sn  Log  pn/Sn 


to  Sn  and  Lg  ratios  increases  with  frequency  (e.g.,  Blandford,  1981;  Bennett  et  al,  1989). 
However,  it  should  be  noted  that  for  the  western  Norway  events,  this  result  is  due  to  the  Pn 
amplitude  being  much  increased  relative  to  the  Lg  amplitude. 

It  is  surprising  that  the  amplitude  ratios  for  the  confirmed  blasts  in  the  same  blasting  site 
have  greater  variance  in  the  amplitude  ratios  than  the  suspected  earthquakes  in  Region  1  and  4 
which  are  distributed  over  a  much  larger  area.  The  blast  variations  may  be  related  to  differences  in 
the  shooting  practice  from  event  to  event,  such  as  in  the  number  and  delay  times  of  the  ripple-fired 
shots.  As  has  been  shown  in  earlier  studies  (Baumgardt  and  Ziegler,  1988)  and  this  study,  the 
complexity  of  the  spectral  modulations  vanes  significantly  from  shot  to  shot,  although  these 
vanations  should  be  the  same  for  all  phases.  Sometimes,  the  nulls  introduced  by  the  modulations 
are  very  deep  and  fall  below  the  noise  level,  particularly  for  the  Pn  phases  which  for  some  events 
have  lower  signal-to-noise  ratios  in  the  8  to  16  Hz  band  than  the  Sn  and  Lg  waves.  Notice,  for 
example,  that  the  F.X9  spectrum  in  Figure  8b  has  a  very  deep  null  at  about  8  Hz  in  the  Pn 
spectrum,  which  is  the  first  spectrum  above  the  noise  spectrum  at  the  bottom.  This  null  is  less 
deep  in  the  case  of  the  l.g  spectrum  at  the  top.  Perhaps  the  Pn  amplitudes  are  more  affected  by 
variations  in  the  number  and  spacing  of  the  deep  nulls  that  result  in  more  noise  contamination  than 
the  Sn  and  l.g  amplitudes,  causing  the  observed  variation  in  amplitude  ratios.  No  such  variation  is 
seen  for  the  presumed  underwater  explosions  in  Region  5  because  they  have  more  similar  spectral 
modulations  than  the  confirmed  blasts  on  land,  and  the  signal-to-noise  ratios  of  Pn  are  high  enough 
that  the  nulls  do  not  drop  into  the  noise. 

Another  possibility  is  that  differences  in  the  time  and  space  distribution  of  the  shots  in  the 
ripple  fire  pattern  may  actually  produce  varying  amounts  of  shear-wave  energy  compared  with 
eompressional-wave  energy.  However,  examination  of  the  points  on  the  scatter  plots  in  Figure  13 
for  the  confirmed  blasts  shows  that  there  is  more  correlation  between  the  amplitude  ratios  with  the 
Pn  amplitudes  than  with  the  Sn  or  l.g  amplitudes.  This  is  not  consistent  with  the  amplitude-ratio 
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variations  being  caused  by  the  variations  in  the  amount  oi  Sn  and  Lg  wave  energy  produced  by  the 
different  blasts. 

In  summary,  we  find  that  regional  Pn'Sn  and  Pn  Lg  ratios  discriminate  well  between 
known  blasts,  suspected  blasts,  and  known  earthquakes  in  western  Norway,  and  that  this 
discriminant  does  not  seem  to  be  significantly  related  to  regional  propagation  effects.  We  observe 
very  consistent  ratio  values  for  our  earthquake  population,  and  the  variations  observed  in  the  blast 
population  seem  to  be  caused  by  ripple-fire  effects.  We  also  find  that  the  observed  differences  in 
the  ratios  correlate  more  with  variations  in  the  Pn  amplitudes  than  with  the  Sn  and  Lg  amplitudes, 
and  that  our  populations  of  blasts  and  earthquakes  in  the  2.0  to  2.5  local  magnitude  range  generate 
very  similar  amounts  of  shear-wave  energy.  However,  our  sample  size  is  too  small  to  determine  if 
this  is  a  result  of  selection.  In  trying  to  get  events  of  comparable  magnitude,  determined  by  coda 
(probably  Lg  coda)  duration  magnitudes,  we  may  have  purposely  selected  events  of  comparable  Lg 
excitation.  If  we  nad  selected  events  on  the  basis  of  some  Pn  magnitude,  we  may  have  found  more 
correlation  of  the  amplitude  ratios  with  the  Sn  and  Lg  amplitudes  rather  than  with  the  Pn  amplitude. 

Comparison  of  Lg  Spectral  Ratios 


We  now-  consider  regional  variations  in  the  ratio  of  low-to-high  frequency  energy  in  the  Lg 
spectrum,  a  feature  which  has  been  found  to  effectively  separate  nuclear  explosions  and 
earthquakes  in  the  western  I'nited  States  (Murphy  and  Bennett,  1982:  Bennett  and  Murphy,  1986; 
Taylor  et  al.  1988).  As  mentioned  above,  the  Lg  spectra  have  all  been  corrected  for  instrument  and 
Pn  background  noise  lor  each  of  the  Lg  spectra,  the  ratios  of  the  rms  levels  in  the  2  to  6  Hz  to 


the  6  to  10  Hz  bands  were  computed  as  follows: 


(2  -  6  Hz  ) 
(6  10  U:  ) 


f  igure  M  shows  the  Lg  spectral  ratio  values  by  region.  Previous  discrimination  studies  in 
the  w  estern  Ini  ted  States  found  that  earthquakes  had  low  er  spectral  ratios  than  explosions. 
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indicating  that  earthquakes  have  more  high  frequency  content  than  explosions.  However,  our 
populations  of  earthquakes  and  blasts  almost  completely  overlap  each  other  in  their  spectral  fatios. 
The  presumed  underwater  explosions  in  Region  5  have  somewhat  higher  spectral  ratios  than  the 
nearby  TITA  blasts.  However,  these  events  clearly  have  lower  spectral  ratios  than  the  suspected 
earthquakes  in  Region  1.  Since  the  TITA  blasts  are  in  Region  5,  this  difference  must  be  an 
indication  of  source  differences,  with  the  off-shore  underwater  blasts  having  more  low  frequency 
content  than  the  on-shore  TITA  blasts.  Thus,  these  results  show  that  spectral  ratio  does  not 
discriminate  between  blasts  and  earthquakes,  but  may  serve  a  means  of  distinguishing  underwater 
blasts  from  blasts  on  land. 

As  in  the  case  of  the  amplitude  ratios,  there  seems  to  be  no  significant  regional  dependence 
of  the  Lg  spectral  ratio.  The  earthquakes  in  Region  4  have  similar  spectral  ratios  to  those  in  Region 
1.  Moreover,  ripple-fire  effects  do  not  seem  to  greatly  affect  spectral-ratio  measurements,  since  the 
variance  of  the  spectral  ratios  for  the  TITA  and  BLA  blasts  is  comparable  to  that  of  the  earthquake 
group. 

2.4  DISCUSSION 

The  consideration  of  case-based  methods  has  been  motivated  by  the  failure  of  previous 
seismic-discrimination  research  to  develop  a  set  of  consistent  regional  seismic  discriminants  and  a 
model  which  explains  how  intrinsic  seismic-source  differences  affects  seismic- waveform  features. 
Lacking  such  a  model,  we  would  identify  seismic  events  by  simply  comparing  them  with  historic 
events  of  of  known  identity  and  not  worry  about  the  explanation  of  why  the  waveform  features  of 
explosions  and  earthquakes  differ. 

However,  this  study  has  pointed  out  some  of  the  problems  that  would  be  encountered  in  an 
operational  setting  in  trying  to  use  case-based  methods  to  identify  seismic  events  to  monitor  a  low- 
yield,  nuclear  testban  treaty  using  the  regional  arrays.  Ideally,  a  seismic-discrimination  experiment 
requires  there  to  lie  large  populations  of  known  earthquakes,  economic  blasts,  and  low-yield 
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nuclear  explosions  in  the  same  region  so  that  useful  waveform  features  sensitive  to  source 
differences,  not  propagation-path  differences,  can  be  identified  as  discriminants.  However,  in 
reality,  for  any  given  region,  this  situation  will  almost  never  be  realized.  At  low  magnitudes,  there 
will  be  many  events  detected  by  the  arrays,  both  natural  and  man-made,  which  will  usually  not 
occur  in  the  same  region.  Moreover,  in  any  future  testban  treaty  scenario,  there  will  almost 
certainly  be  no  small  nuclear-explosion  tests  available  for  comparison. 

Our  study  has  shown  how  case-based  methods  may  be  useful,  in  spite  of  all  the 
difficulties.  The  approach  would  involve  first  establishing  a  baseline  of  "normar  seismic  activity 
in  a  region,  characterized  by  the  location  of  historic  events  and  the  nature  of  the  waveform  features 
for  the  events.  As  much  as  possible,  demographic  information  would  be  exploited  (e.g.,  locations 
of  events  relative  to  high-population  centers,  mines,  and  previous  test  sites,  felt  reports,  time  of 
occurrence  of  events)  although  not  relied  on  exclusively.  It  should  be  possible  to  have  a 
population  of  confirmed  explosions  and  other  events  which  are  most  likely  earthquakes,  although  it 
will  be  hard  to  positively  identify  earthquakes  in  areas  of  low  natural  seismicity.  We  have  found  in 
our  study  that  many  of  the  discriminants  appear  to  work  in  identifying  obvious  event  types,  and 
that  the  propagation  effects  in  western  Norway  seem  to  be  essentially  homogeneous  throughout  the 
region.  Earthquakes  can  be  identified  on  the  basis  of  low  P-to-S  ratios  and  lack  of  spectral 
modulations.  Mine  and  off-shore  explosions  usually  have  time-independent  spectral  modulations 
and  high  P-to-S  ratios.  Obvious  event  types  can  be  quickly  identified,  based  on  these  features. 

An  important  question  about  discrimination  raised  by  this  study  is  why  the  spectral  ratio 
discriminant  fails  to  work  as  well  for  separating  earthquakes  and  economic  explosions  in  western 
Norway  as  it  did  in  separating  small  nuclear  explosions  and  earthquakes  in  western  United  States 
(Murphy  and  Bennett,  1986;  Bennett  and  Murphy,  1986),  whereas  the  P-to-S  ratio  seems  to  work 
better  in  western  Norway  than  it  did  in  western  United  States.  Bennett  et  al  (1989)  also  found  that 
P-to-S  ratios  worked  better  in  separating  earthquakes  and  explosions  in  the  Soviet  Union  than  Lg 
spectral  ratio,  although  their  earthquakes  and  explosions  were  not  located  in  the  same  regions.  In 
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our  study,  propagation-paths  effects  do  not  seem  to  have  a  strong  effect  on  the  Lg  spectral  ratios. 
Moreover,  large  variations  in  spectral  ratios  were  observed  for  mine  explosions  in  the  same 
location,  presumably  caused  by  the  ripple-fire  effects.  However,  we  also  found  that  even  for  the 
closely  clustered  earthquakes  in  Region  1,  whose  P-to-S  ratios  were  very  similar,  the  Lg  spectral 
ratios  exhibited  high  variance. 

Lilwall  (1988)  has  offered  an  explanation  of  the  spectral  ratio  discriminant  in  terms  of  the 
effects  of  evanescent  waves  or  S  on  Lg-wave  generation.  Such  waves  would  be  be  more 
strongly  generated  by  shallow  explosions  near  the  surface  than  deeper  earthquakes,  and  thus,  the 
S  waves  would  build  up  the  low  frequency  levels  in  the  Lg  waves  from  explosions  relative  to 
earthquakes.  Such  effects  do  not  seem  to  be  present  in  our  blast  group.  In  fact,  the  underwater 
explosions  in  Region  5  seem  to  have  slightly  lower  values  of  R  than  other  events  (including  mine 
explosions),  indicating  that  the  higher  frequencies  in  Lg  are  enhanced  for  these  events. 

Taylor  et  al  (1988)  suggests  that  the  presence  of  a  low-Q  region  in  the  upper  crust  of  the 
western  United  States  might  explain  why  explosions,  which  are  shallow  and  occur  in  the  low  Q 
region,  have  less  high  frequencies  than  earthquakes,  which  might  occur  below  the  highly 
attenuative  region.  "  rhaps  no  such  low-Q  region  exists  in  western  Norway  thus  explaining  why 
the  Lg  spectra  of  both  earthquakes  and  explosions  are  similar.  The  performance  of  the  P-to-S  ratio 
may  be  more  a  result  of  intrinsic  source  differences,  i.e.,  explosions  have  less  shear  energy 
compared  to  compressional  energy  at  high  frequency  than  earthquakes.  Our  P-to-S  ratio  data 
seems  to  correlate  more  with  the  P-wave  levels  than  with  the  S-wave  levels.  This  suggests  that 
explosions  and  earthquakes  generate  the  same  amount  of  shear-wave  energy,  but  that  explosions 
produce  more  compressional-wave  energy  than  earthquakes.  However,  this  may  be  more  a  result 
of  event  sampling  based  on  the  Bergen  local  magnitude,  as  was  discussed  above.  Differences  in 
depth  of  focus  may  also  be  partly  responsible.  The  mine  explosions  apparently  occurred  at  the 
surface  (Svcin  Mykkeltveit,  personal  communication)  whereas  the  earthquakes  were  deeper. 
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Perhaps  near-surface  effects,  such  as  PS  conversions  and  S*  generation,  may  build  up  the  shear 
wave  energy  for  the  earthquakes. 

This  study  has  shown  how  seismic  events  can  be  identified,  relative  to  case  events,  by 
means  of  a  systematic  assessment  of  the  similarities  and  differences  between  waveform 
characteristics,  alwavs  keeping  in  mind  the  possible  effects  of  propagation-path  differences.  We 
have  found  that  array-averaged  spectra  and  incoherent  beams,  computed  using  regional-array  data, 
provide  robust  estimates  of  waveform  features  important  for  seismic  discrimination.  We  have 
shown  that  the  explosions  and  earthquakes  in  a  small  region  of  western  Norway,  recorded  at  the 
NORESS,  can  be  well  separated  on  the  basis  of  high-frequency  amplitude-ratio  discriminants,  but 
that  the  Tg  spectral  ratio  discriminants  do  not  separate  as  well.  This  seems  to  agree  in  general  with 
other  regional  discrimination  studies  in  Scandinavia  and  Eurasia. 

This  study  has  also  demonstrated  that  the  ripple-fire  discriminant  based  on  the  observation 
of  time-independent  modulations,  proposed  by  Baumgardt  and  Ziegler  (1988)  and  Hedlin  et  al 
1 1989).  may  prove  to  he  highly  useful  in  identifying  many  unknown  regional  events  which  are 
produced  by  unreported  economic  blasting.  Apparently,  at  the  low-magnitude  monitoring  level 
required  for  a  low  yield  or  comprehensive  testban  treaty,  many  such  events  will  be  detected  by  the 
regional  arrays. 

Because  discriminants  seem  to  be  inconsistent  in  their  performance  and  that  we  are  not  sure 
how  u>  correct  discriminants  for  propagation-path  effects,  we  have  advocated  the  use  of  a  case- 
based  approach  to  regional  seismic  discrimination.  In  this  approach,  events  should  be 
characterized  in  a  step-by-step  fashion,  testing  individual  discriminants,  and  always  trying  to  relate 
the  signal  characteristics  of  unidentified  events  to  events  of  known  identity  in  the  same  or  similar 
regions  Because  of  me  curicm  lack  of  a  unifying  theory  for  regional  seismic-event  identification, 
the  case  based  approach  may  be  the  only  wav  to  reliably  characterize  and  identify  seismic  events. 
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SKCTION  3.0 


RIITLK-HRK  DKTKCTION 


3.1  INTRODl  C'TION 

In  the  study  described  in  the  previous  section,  we  have  found  that  many  events  recorded  at 
NORKSS,  w  hich  had  prev  iously  been  suspected  to  he  earthquakes  because  they  were  not  reported 
to  be  mine  blasts,  had  spectra  which  revealed  time-independent  spectral  modulations  or 
scalloping."  Several  previous  studies  (Bell  and  Alexander.  1977;  Baumgardt  and  Ziegler.  19XX; 
Smith.  19X9.  Hedlin  et  al,  19X9.  1990)  have  shown  that  this  is  the  distinctive  signature  of  multiple 
seismic  events,  or  "ripple  firing,"  which  is  a  very  common  practice  in  most  commercial  blasting 
activity  .  We  concluded  that  the  events  were  more  likely  to  be  ripple-fired  blasts  than  earthquakes, 
since  the  spectral  scalloping  was  so  strongly  pronounced  and  other  discriminants  seemed  to  be 
more  consistent  with  the  events  being  explosions.  For  this  reason,  we  believe  that  most  regional 
blasting  activity  can  be  identified  if  this  ripple-fire  pattern  can  be  observed. 

The  ISKIS  system  design  will  include  a  subsystem  which  will  automatically  identify  time- 
independent  spectral  modulations  in  regional  seismic  data  and  will  serve  as  a  primary  discriminant 
for  the  identification  of  blasting  activity.  Two  methods  have  been  suggested  in  the  literature  for 
automatically  detecting  spectral  scalloping.  Baumgardt  and  Ziegler  (198X)  showed  that  time- 
independent  spectral  scalloping  produces  quefrency-independent"  peaks  in  power  cepstra,  which 
are  Fourier  transforms  of  the  modulated  spectra.  As  part  of  the  development  of  the  IAS.  a  svstem 
called  the  Multiple  Hvcnt  Recognition  System  (MHRSY)  was  developed  which  automatically 
detected  quefrenev  independent  cepstral  peaks  in  cepstra  computed  for  two  or  more  regional 
phases  associated  with  an  event.  This  work  was  described  in  an  unpublished  report  by  Baumgardt 
and  Ziegler  (  19X9)  The  second  method,  described  by  Hedlin  et  al  (1990),  also  uses  cepstra 
computed  for  entire  reition.il  coda  derived  from  the  Fourier  transform  of  the  time  frequence 
spectrograms  or  sonograms  of  the  c<xla  However,  their  method  looks  for  high  cepstral  levels. 
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relative  to  a  threshold  established  from  unmodulated  spectra  of  earthquakes,  as  the  criterion  for  the 
identification  of  modulated  spectra  due  to  ripple-fire.  1’heir  method  is  similar  to  using  cepstral 
variance  as  the  criterion,  as  suggested  as  a  secondary  criterion  in  the  MERSY  approach 
(Baumgardt  and  Ziegler.  19X9).  High  cepstral  values,  due  to  spectral  modulations,  would  also 
produce  high  cepstral  variance. 

In  (Mir  view  ,  it  is  essential  that  the  time  independence  of  spectral  modulations  be  identified, 
since  other  effects  can  produce  high  cepstral  value  or  high  cepstral  variance,  l  or  example,  in  the 
Pn  coda,  there  may  be  crustal  reverberations  and  mode  conversions  from  high  impedance-contrast 
interfaces  in  the  crust  beneath  the  source  and/or  receiver  which  could  produce  complex  spectral 
modulations  and  high  cepstral  levels.  Similarly,  Sn  coda  may  also  contain  shear-wave 
reverberations  from  the  same  crustal  interfaces  which  could  produce  complex  modulations  in  the 
.S'/i-coda  spectra  and  would  also  produce  high  cepstral  levels.  Thus,  the  presence  of  high  cepstral 
levels  in  both  the  Pn  and  Sn  codas  does  not  by  itself  necessarily  imply  that  the  event  is  ripple-fired. 
Only  the  observation  of  the  same  modulation  pattern  throughout  the  entire  coda  can  unambiguously 
be  used  to  identify  ripple  tired  explosions. 

Bennett  et  a!  ( I990i  have  criticized  the  cepstral  method  by  pointing  out  that  cepstral  peaks 
can  be  observ  ed  in  the  P  codas  of  earthquakes  in  western  US.  However,  it  is  important  to  realize 
that  it  is  the  quefrency  independence  of  the  cepstral  peaks,  not  the  presence  of  the  cepstral  peaks 
themselves,  which  is  key  to  identifying  ripple-fired  events.  Stump  and  Reinke  (19XX)  have 
pointed  out  that  spectral  modulations  may  be  less  clear  if  there  are  significant  errors  in  the  delay 
times  of  the  blasting  timers  used  in  a  delayed  shooting  sequence.  Moreover,  multiple  shots  set  off 
in  complex  spatial  patterns  would  produce  complex  modulations  in  the  spectra  (e.g..  Smith.  19X9). 
However,  whatever  modulations  arc  produced,  even  if  not  strong,  will  still  be  time-independent  in 
the  ease  ot  ripple  filing  t  hus,  even  weak  cepstral  peaks  produced  by  weak  spectral  modulations, 
which  might  tail  the  tests  suggested  by  lledlin  el  al  (  19X9).  would  still  be  quefrency  independent 
and  could  be  a  strong  indication  ot  ripple  firing,  if  they  can  be  detected  reliably . 
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f  or  the  ISHIS  system,  we  chose  to  adapt  the  MHRSY  approach  from  the  IAS,  since  the 
method  has  proven  to  he  effective  on  a  large  database  of  events  used  to  test  the  IAS.  The  method 
has  been  discussed  in  detail  by  Baumgardl  and  Ziegler  ( 19X9i  with  numerous  examples  presented. 
However,  since  that  report  has  not  yet  been  generally  published,  a  description  of  the  method  and 
some  examples  are  presented  below 

3.2  RIPPI.K-FIRK  DKTKCTION  METHOD 

The  Multiple  Event  Recognition  System  (MHRSY)  was  designed  to  be  part  of  the  post¬ 
detection,  event  characterization  processing  of  the  Intelligent  Array  System  (IAS).  The  design 
requirements  for  this  system  are  described  in  the  MHRSY  requirements  document  (ENSCO, 
19X8).  Essentially,  MHRSY  automates  the  method  suggested  by  Baumgardt  and  Ziegler  (1988)  of 
identifying  spectral  modulations  by  computing  Fourier  spectra  on  each  phase  associated  with  an 
event,  then  computing  the  cepstra  by  Fourier  transforming  the  spectra,  and  then  comparing  the 
resulting  cepstra  to  find  two  or  more  peaks  which  have  the  same  quefrency. 

Fourier  Cepstral  Analysis 

The  methodology  for  cepstral  analysis  has  been  described  in  detail  by  Baumgardt  and 
Ziegler  (1988).  The  procedure  simply  involves  taking  the  Fourier  transform  of  the  base-ten 
logarithm  of  the  spectral  density,  computed  by  Fourier  transforming  the  windowed  time-series  of 
the  detected  phases.  These  spectra  were  computed  by  averaging  the  spectra  computed  for  each 
channel  of  the  array,  following  the  method  described  by  Baumgardt  and  Ziegler  (1987).  The 
cepstrum  is  closely  related  to  the  auto  correlation  function,  which  is  also  the  Fourier  transform  of 
the  power  spectrum.  The  mam  difference  between  the  cepstrum  and  the  auto-correlation  function 
is  that  the  logarithm  of  the  spectra!  density  is  taken  in  computing  the  cepstrum.  which  effectively 
whitens  the  spectrum  thus  producing  sharper  cepstral  peaks. 
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Before  computing  the  second  Fourier  transform,  any  systematic  trends  must  be  removed 
from  the  spectra  for  if  not,  the  trend  itself  will  produce  spurious  cepstral  peaks.  Originally,  the 
spectra  were  assumed  to  decay  linearly  with  frequency.  A  linear  trend  was  fit  by  least-squares  to 
the  spectra  above  2  Hz,  and  then  the  linear  trend  was  subtracted  from  the  data.  However,  we 
found  that  most  of  the  spectra  had  non-linear  trends,  as  shown  in  the  example  in  Figure  15.  The 
bottom  show  s  the  spectrum  of  an  /.g  wave  from  a  ripple-fired  explosion  fit  by  least  squares  with  a 
quadratic  curve.  The  top  plot  shows  the  same  spectrum  with  the  quadratic  trend  removed,  which 
has  essentially  a  fiat  slope  on  which  the  spectral  modulation  can  be  seen  clearly.  Thus,  the  linear 
trend  removal  has  been  replaced  with  a  quadratic  trend  removal,  a.>  shown  in  Figure  15. 

Maximum  Entropy  Cepstra 

One  of  the  well-known  problems  associated  with  cepstrum  analysis  is  that  the  sharpness  of 
cepstral  peaks,  defined  by  the  maximum  [tower  of  the  peak  and  its  handwidth,  is  controlled  by  the 
bandwidth  of  the  spectrum.  Also,  the  shorter  the  time  delays  of  delayed  explosions,  the  less  the 
number  of  cycles  in  the  modulation  within  the  bandwidth  of  the  spectrum.  For  example,  the 
sampling  rate  of  high  frequency  NORFiSS  data  is  40  Hz  which  gives  a  20  Hz  N'yquist.  Thus,  the 
shortest  delay  which  can  resolved  with  this  data  is  25  ms.  which  would  give  a  modulation  of 
exactly  20  1  lz.  In  order  to  resolve  shorter  delays,  it  would  be  necessary  to  use  broader-band  data. 

Another  problem  is  that  even  for  delays  greater  than  50  ms.  conventional  Fourier  cepstrum 
analysis  involves  truncating  the  spectrum  at  some  high  frequency  cut-off  point,  usually  the 
Nvquist  frequency.  This  could  cut-off  a  spectral  modulation  in  midcycle  and  produce  a  spurious 
peak  at  a  diflercnt  quefrency  than  the  actual  quefrenev  of  the  modulation.  Moreover,  aliasing  and 
noise  effects  at  frequencies  near  the  Nvquist  frequency  can  produce  other  spurious  peaks. 


44 


in 


Spectral  Detrend 


r 


FKilRE  15:  Kxample  of  the  fit  of  a  second-order  polynomial  to  the  spectrum 
recorded  at  NORKSS.  The  dashed  line  indicates  the  pre-Pn  noise  spectrum.  The 
top  plot  is  the  spectrum  after  the  quadriatic  trend  has  been  remo\ed. 


We  implemented  and  experimented  with  a  maximum  entropy  method  (MEM)  algorithm  for 
computing  the  cepstrum  in  MHRSY  in  order  to  eliminate  cepstral  peaks  due  to  noise  and  processing 
anifact.  We  were  motivated  by  the  fact  that  the  MEM,  in  spectral  estimation,  produces  sharper 
spectral  peaks  because  it.  in  effect,  extrapolates  the  auto  correlation  function  to  time  lags  greater 
than  the  time- window  length.  In  our  use  of  MEM  for  cepstrum  analysis,  we  treated  the  log- 
amplitude  spectrum  as  if  it  were  a  time-series  and  computed  the  "MEM  cepstrum"  which  would 
have  the  effect  of  extrapolating  the  frequency-band  window.  Thus,  we  wo  Id  expect  that  MEM 
cepstrum  analysis  would  broaden  the  bandwidth,  extrapolate  spectral  modulations  to  frequencies 
beyond  Nyquist.  thus  eliminating  spurious  peaks  due  to  truncation  and  sharpen  the  peaks  caused 
by  the  spectral  modulation. 

T  he  Burg  (1967)  maximum  entropy  power  spectrum  formalization  was  adopted  and  the 
MEM  process  was  implemented  using  the  algorithm  described  by  Anderson  (1974),  except  that  the 
process  was  applied  to  the  spectrum  instead  of  to  the  time-series.  The  MEM  power  cepstrum  is 
estimated  by 
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where  the  quefrency.  t.  is  limited  to  the  interval,  defined  by  -1/(21) f)  </<  1/(2D/)  where  D/is  the 
frequency  sampling  interval  and  amn  are  determined  by  the  equation. 
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This  equation  gives  the  elements  of  a  predictive  filter,  (I ,-ami . i*mn  )•  derived  from  the 

spectral  auto  correlation  function  with  frequency  lag  /.  or  /  /.  Pm  is  the  output  power  of  the  m  +  l 
long  predictive  filter,  given  above,  f  or  m-().  Pa  is  estimated  by 
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The  MEM  algorithm  iteratively  solves  the  above  matrix  equation,  using  the  algorithm  of 
Anderson  (1974),  except  that /should  be  substituted  for  t. 

Spectral/Cepstral  Statistics 

MERSY  also  computes  the  following  statistically  based  features  on  the  spectra  and  cepstra. 

The  variance  estimate  is  a  measure  of  the  spectral  or  cepstral  "width"  or  "variability"  around 
the  mean  value.  Let  xj  represent  the  jth  spectral  or  cepstral  value  and  x  the  mean  value,  then  the 
variance  is  defined  as 
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where  N  is  the  number  of  points  in  the  spectrum  or  cepstrum.  The  variance  is  also  known  as  the 
second  moment. 


The  skewness  characterizes  the  degree  of  asymmetry  of  a  function  or  distribution  about  its 
mean  value.  It  is  represented  as  the  third  moment  expressed  as 
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where  s  -  six /  ...  xn)  is  the  standard  deviation  of  the  spectral  or  cepstral  distribution.  An  estimate 
of  the  standard  deviation  is 


...xN  )  =  JVar(xx...xN  ) 


The  kuriosis  gives  a  measure  of  the  "peakedness"  or  "flatness"  of  the  spectrum  or  cepstrum 
relative  to  a  normal  distribution  function.  This  feature,  also  known  as  the  fourth  moment,  is 
expressed  as 


Kurt  u, . in  )  = 
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The  -  3  term  makes  the  value  of  Kurt  zero  when  the  spectrum  or  cepstrum  is  normally 
distributed. 

These  features  may  serve  two  purposes.  First,  they  measure  the  shape  of  the  spectra  and 
cepstra  which  may  be  useful  for  discrimination.  For  example,  Pulli  and  Dysart  (1987)  showed  that 
the  variance  of  the  cepstrum  discriminated  between  several  earthquakes  and  mine  explosions  in 
Scandinavia  recorded  by  NORKSS.  Apparently,  the  mine  blasts  had  higher  variance  than  the 
earthquakes  as  a  result  the  spectral  modulations  in  the  former  caused  by  the  ripple-fire  effects. 
Buumgardt  and  Ziegler  ( 19X9)  also  found  that  populations  of  earthquakes  and  explosions  could  be 
separated  on  the  basis  of  variance,  skew,  and  kurtosis  of  the  cepstra  for  region,  phases.  Thus, 
these  statistics  can  serve  as  secondary'  features  for  identification  of  ripple-fire  in  addition  to  the 
observ  ation  of  cepstral  peaks. 

The  spectral  variance  is  also  useful  for  determining  the  number  m  of  coefficients  utnn 
which  must  be  included  in  the  maximum  entropy  cepstrum  predictive  filter  for  n  frequency  points. 
The  number  of  coefficients  depends  on  the  complexity  of  the  cepstrum.  which  is  partly  measured 


4X 


by  the  logarithm  of  the  Fourier  cepstral  variance,  Icvar.  The  following  empirically  derived  rules 
are  used  to  set  the  number  of  coefficients 

Icvar  -4. 1  <  -4.4  — >  m=5() 

Icvar  >  -4.4  and  Icvar  <  -4.1  — >  m  =  70 
Icvar  >-4.1  — >  m  =  80. 

In  addition  to  these  statistics,  an  estimate  of  the  spectral  signal-to-noise  ratio  (ssnr)  is  also 
made.  The  spectral  snr  is  defined  as  the  bandlimited  ratio  of  the  rms  signal  to  the  rms  noise,  or 
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where  i  I  and  Mfare  the  first  and  last  frequency  points  of  the  band,  5/  is  the  spectrum  of  the 
detected  ph.  se  and  Nm  is  the  the  background  noise  sample  to  the  first  of  a  detection  group.  The 
summation  is  taken  over  the  frequency  band  from  1.88  to  19.7  Hz. 

Single-Array  Peak  Analysis 

Peaks  are  identified  in  the  Fourier  and  MEM  cepstra  and  analyzed  to  determine  if  they  are 
caused  by  spectral  modulations  caused  by  ripple-fire.  Using  the  method  suggested  by  Baumgardt 
and  Ziegler  ( 1988),  multiple  explosions  are  identified  based  on  the  appearance  of  a  consistent  set  of 
significant  peaks  having  the  same  quefrency  in  the  cepstra  of  all  the  phases  associated  with  the 
event.  Significant  peaks  are  defined  as  cepstral  peaks  whose  amplitudes  exceed  some  threshold. 
Consistent  peaks  are  those  which  appear  in  two  or  more  phases,  but  not  in  the  noise  cepstra.  The 
specific  criteria  for  the  peaks  to  be  significant  and  consistent  are  described  in  ENSCO  (1989).  In 
brief,  two  of  the  following  three  conditions  must  hold  for  the  event  to  be  classified  as  a  multiple 
event: 
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(1)  One  or  more  consistent  significant  peaks  appearing  in  the  Fourier  cepstra  of  two  or 
more  phases  associated  with  an  event,  hut  observed  in  the  noise  cepstra. 

(2)  Same  as  above  for  MEM  cepstra. 

(3)  The  cepstral  variance  must  exceed  a  minimum  threshold. 

Multi-Array  Analysis 

If  two  or  more  arrays  detect  a  mine  explosion  with  ripple-fire,  there  should  be  consistency 
between  the  spectral  modulations  observed  at  all  the  arrays  as  well  as  between  phases  detected  at 
one  array.  In  order  to  detect  consistent  modulations  at  more  that  one  array,  stack  cepstra  and 
spectra  are  computed  for  each  array  by  averaging  the  cepstra  of  all  the  phases  associated  with  the 
event.  Thus,  there  is  one  stack  spectrum  and  cepstrum  for  each  array.  (Note:  Stack  cepstra  are 
computed  by  averaging  together  the  cepstra  for  each  phase;  they  are  not  computed  from  the  stack 
spectra.  Stack  spectra  are  computed  for  the  purpose  of  computing  signal -to- noise  ratios  and  other 
spectral  features. ) 

The  stack  cepstra  for  each  array  are  treated  in  the  same  manner  as  the  individual  phase 
cepstra  in  single  array  processing.  Significant  peaks  are  identified  on  each  of  the  stack  cepstra  for 
each  array  w  hich  have  the  same  quefrency  at  two  or  more  arrays. 

3.3  MERSY  PROCESSING  EXAMPLES 

In  this  section,  we  show  some  examples  of  MERSY  processing  results  for  events 
processed  in  the  IAS.  Figure  16  shows  an  example  of  a  mine  explosion  in  Estonia  processed  by 
the  MERSY  system.  The  top  left  window  shows  three  filtered  beam  traces  with  the  identified 
phases  Pn.  P  S n.  and  l-K  marked  The  right  three  windows  show  the  spectra,  at  top.  the  Fourier 
l EFT)  cepstra.  in  the  middle,  and  the  MEM  cepstra.  at  the  bottom.  Each  of  the  spectra  were 
computed  in  six  second  w  indows,  beginning  at  the  marked  phase  times,  and  have  been  corrected 
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ru.ijut  10.  t.xampie  01  a  vilkjvy  display  lor  a  presumed  Estonia  mine  blast  using  the  X  window  graphics. 
Version  10.4.  The  upper  left  plot  show's  three  beam  waveforms  (vertical,  horizontal,  incoherent).  The  left  three 
panels  show  the  spectra,  FFT  cepstra,  and  MEM  cepstra  for  the  four  phases  detected  and  associated  with  the 
event  and  the  pre-Pn  noise.  The  vertical  lines  on  the  cepstra  indicate  the  peaks  picked  bv  the  MERSY  algorithm. 
The  lower  right  window  is  a  textual  display  of  the  processing  results. 


for  instrument.  The  spectra  are  vertical  component  spectra  averaged  across  the  array.  The  array 
data  is  sampled  at  40  Hz,  giving  a  Nvquist  frequency  of  20  Hz.  The  spectra  and  the  cepstra  are 
displayed  for  each  phase  associated  with  the  event,  along  with  the  pr eTn  noise,  labeled  as  NS. 
Note  that  the  spectra  and  cepstra  have  been  displaced  upward,  relative  to  the  Pn  spectrum  and 
cepstrum,  for  display  purposes.  1'he  window  on  the  lower  left  gives  a  textual  description  of  the 
results  of  the  MERS  Y  processing,  including  a  statement  of  the  number  of  peaks  identified  in  the 
FFT  and  MEM  cepstra,  the  average  cepstral  variance,  and  the  low  frequency  cutoff  of  the  band. 
The  table  gives  the  values  of  the  spectral  signal-to-noise  ratios,  the  spectral  and  cepstral  variances, 
and  the  number  of  prediction  error  coefficients  used  in  computing  the  MEM  cepstra  for  each  phase. 
Finally,  the  actual  quefrencies  of  the  consistent  peaks  are  written,  which  correspond  to  the 
dominant  delay  time  of  the  multiple  events. 

For  the  case  of  the  Estonia  mine  explosion  in  Figure  16,  a  clear  modulation  pattern  can  be 
seen  in  the  spectra  for  all  the  phases.  The  nulls  in  the  spectral  scalloping  are  spaced  apart  by  about 
0.25  Hz,  which  corresponds  to  a  delay  time  of  about  400  ms.  On  the  cepstral  plots,  vertical  lines 
have  been  plotted  through  the  peaks  which  were  identified  by  MliRSY  as  being  significant  and 
consistent.  MI  ES'!  picked  the  primary  peak  on  the  Fourier  (FFT)  cepstra  at  0.4  seconds,  as 
expected,  in  addition  to  a  secondary  peak  at  about  IX  seconds.  Note  that  the  spectral  modulations 
are  not  observed  in  the  noise  spectra  nor  do  cither  of  the  cepstra!  peaks  appear  in  the  noise  cepstra. 

The  MFM  cepstrum  in  the  bottom  right  shows  two  consistent  peaks,  which  were  identified 
by  MI  ES')  and  indicated  by  the  vertical  lines.  These  peaks  are  significantly  sharper  than  those  in 
the  Ilf  cepstra  and  the  variance  of  the  rest  of  the  cepstrum  has  been  much  reduced.  However,  the 
quefrencies  of  these  two  peaks  (  150  and  425  ms)  are  shifted  relative  to  the  quefrencies  (175  and 
400  ms)  of  the  MFM  cepstral  peaks.  Also,  note  that  the  MEM  peaks  are  not  as  consistent  as  the 
FI  1  peaks 
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Peak  shitting  of  this  kind  is  a  well  known  problem  in  maximum  entropy  power  spectral 
estimation  studies  (e  g.,  ('hen  and  Stegen,  1974:  Kaveh  and  Eipper.  19X3).  The  causes  of  such 
peak  shifting  include  ( 1  >  nonstationarity  of  the  input  time  series.  (2)  variations  in  the  initial  phase 
of  the  estimated  time-senes,  and  (3)  noise  contamination.  Baumgardt  and  Ziegler  ( 19X9)  have 
discussed  the  causes  of  cepstral  peak  shifting  and  have  suggested  some  possible  solutions  to  the 
problem. 

One  well-known  MEM  problem  we  have  avoided  in  the  design  of  MKRSY  is  peak 
splitting.  This  problem  arises  when  too  many  terms  are  used  in  the  predictive  filter  (Chen  and 
Stegen,  1974).  As  in  spectral  estimation,  we  have  also  observed  cepstral  peak  splitting  in  MEM 
cepstra  if  too  many  predictive  error  coefficients  are  used  in  the  frequency  domain.  MERSY 
determines  the  number  of  terms  to  use  from  the  variance  of  the  EFT  cepstrum:  i.e.,  the  higher  the 
variance,  the  larger  the  number  of  terms.  The  rule  used  to  determine  the  number  of  coefficients 
from  the  cepstral  variance  was  derived  empirically.  Many  cases  were  run,  under  various  signaJ-to- 
noise  ratios,  in  order  to  determine  the  minimum  number  of  coefficients  to  use  to  provide  high  peak 
resolution,  but  not  too  high  to  cause  cepstral  peak  splitting.  In  all  the  cases  we  have  analyzed 
using  these  rules,  we  have  seen  no  cases  of  peak  splitting. 

Baumgardt  and  Ziegler  ( 19XX)  have  pointed  out  that  .here  are  spurious  peaks  observed  in 
cepstra  due  to  processing  artifact,  particularly  at  quefrencies  near  the  inverse  of  the  Nyquist 
frequency  (0.3  sec).  In  Figure  16.  a  large  negative  peak  (i.e.,  a  tough)  can  be  seen  near  0.5 
seconds,  which  we  believe  is  caused  by  truncation  of  the  spectrum  at  Nyquist.  We  have  observed 
that  the  MEM  cepstrum  «'«ns  to  eliminate  these  low-quefrency  peaks  due  to  processing  artifact. 
The  MERSY  rules  do  not  permit  the  picking  of  peaks  or  troughs  at  this  low  a  quefrency. 

Figure  17  shows  the  steps  involved  in  multi-array  processing.  In  this  case,  data  were 
processed  from  a  mine  explosion,  V 1 2.  near  1  .enmgrad.  ESSR,  recorded  at  both  the  NORESS  and 
FINESS  regional  arrays  The  NORESS  spectra  and  cepstra  are  plotted  on  the  left  and  those  for 


FKfL'RE  17:  An  example  of  two-array  processing  of  a  Leningrad  mine  blast  (V12).  The  left  and  center  three 
windows  are  the  processing  results  for  the  detections  at  NORESS  and  FINKSS.  The  plots  on  the  right  were 
obtained  by  averaging  the  spectra  and  cepstra  of  the  phases  detected  at  each  array. 


FINESS  in  the  middle.  On  the  right,  the  stack  spectra  and  cepstra  are  shown.  In  addition  to  the 
stack  signal  spectra  and  cepstra.  the  noise  spectra  and  cepstra  at  each  array  are  plotted  and  labeled 
as  NS-NOR  and  NS  I  IN  Modulations  are  clearly  observable  in  the  spectra  at  both  arrays  which 
are  identical.  Cepstral  peaks  were  also  identified  in  all  three  FIT  cepstra,  although  the  number  of 
peaks  varies.  In  all  cases,  a  peak  at  about  0. 1 5  seconds  has  been  picked.  The  MEM  cepstra  placed 
the  peaks  at  a  somew  hat  higher  quefrency,  or  about  0.2  seconds.  Also,  there  appears  to  be  more 
variability  in  the  location  of  the  MEM  peaks. 

Figure  IS  shows  another  example  of  a  multi  array  processing  result  for  a  V4  mine 
explosion,  which  is  also  near  Leningrad.  This  case  is  an  example  of  a  long-period  spectra! 
modulation  which  has  been  truncated  at  the  Nyquist  frequency.  One  cycle  appears  to  start  at  near 
dc  and  terminate  at  about  1  ^  Hz,  where  the  first  null  can  be  observed.  Then,  a  second  modulation 
begins.  We  expect  that  this  second  modulation  should  have  the  same  period  as  the  first  and  extend 
from  1 4  to  26  1 1/.  However,  the  conventional  EFT  cepstrum  would  make  a  peak  corresponding  to 
the  modulation  from  1 4  Hz  to  20  Hz.  which  gives  a  delay  time  of  about  0.14  seconds.  In  the  FIT 
stack  cepstra,  a  number  of  peaks  near  this  quefrency  are  apparent,  in  particular  the  F1NESS  stack 
cepstra.  The  MEM  cepstra  only  have  one  peak  at  both  NORESS  and  F1NESS  at  about  0.1 
seconds,  which  is  close  to  the  correct  value.  Evidently,  the  MEM  cepstrum  appears  not  to  have  the 
spurious  peaks  caused  by  truncation  of  the  spectrum  at  Nyquist. 

Figure  Id  shows  a  MERSY  result  for  NORESS  and  ARCESS  data  for  an  event  which  had 
no  modulations  which  would  produce  cepstral  peaks  with  quefrencies  less  than  one  second. 
(Note:  close  examination  of  these  spectra  suggests  that  there  may  be  small,  high-quefrenev 
modulations,  which  would  be  consistent  with  multiple  events  with  delays  in  excess  of  one  second. 
These  may  be  caused  b\  afiershtvks  in  the  earthquake  coda.  However,  these  large  delay  times  are 
outside  the  range  we  would  expect  tor  typical  mine-blast  ripple-fire  delays.)  This  event,  which 
was  located  in  Norway  between  the  NORESS  and  ARCESS  arrays,  did  not  occur  near  any  known 
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FIGURE  18:  Example  of  NORESS/FINESS  multi-array  processing,  MERSY  display  for  a  presumed  mine 
explosion.  Notice  the  spurious  FFT  cepstral  peaks,  particularly  at  FINESS,  caused  by  the  truncation  of  the 
spectrum  at  Nyquist  in  the  middle  of  a  spectral  modulation. 


J  j  tyrrff.i  for  f  vrnt  | 


57 


FIGURE  19:  Example  of  NORESS/ARCESS  processing  of  a  presumed  earthquake.  The  spectra  are  essentially 
flat  with  no  strong,  consistent  modulations  apparent  with  quefrecies  less  than  1.0  second. 


mine  site  nor  was  it  identified  as  a  mine  blast  in  the  regional  bulletins  of  the  University  of  Bergen. 
We  conclude  that  this  event  is  an  earthquake. 

In  the  previous  section,  we  pointed  out  a  group  of  events  which  occurred  off  the  southern 
coast  of  Norway  which  had  strong  spectral  modulations.  We  argued  that  these  events  may  have 
been  underwater  blasts,  and  that  the  modulations  were  caused  by  source  multiplicity  produced  by  a 
blast-induced  bubble  pulse.  The  locations  of  three  other  examples  of  off-shore  events,  which  we 
have  processed  through  MHRSY  and  may  be  underwater  blasts,  are  shown  on  a  map  in  Figure  20. 
IJWB1  occurred  in  the  Gulf  of  Bothnia,  UWB2  in  the  North  Sea,  and  UWB3  in  the  Baltic  Sea. 
Figure  21  shows  the  spectra  and  eepstra  for  these  events,  computed  from  spectra  of  phases 
detected  at  NORESS.  In  each  case  a  very  clear  spectral  modulation  pattern  can  be  seen  in  all 
phases  which  produces  very  sharp  cepstral  peaks.  As  a  general  rule,  these  events  appear  to  have 
delay  times  on  the  order  of  3(H)  to  600  ms.  These  modulations  are  so  clear  that  they  seem  to  be 
caused  by  the  interference  of  two  pulses  delayed  by  about  0.5  seconds.  Baumgardt  and  Ziegler 
( 19S9)  showed  how  such  modulations  are  consistent  with  explosive  charges  of  about  2500  lbs 
detonated  at  water  depths  of  about  300  ft. 

3.4  CONTI.  UNIONS 

Our  initial  assessment  of  the  MHRSY  system  as  part  of  the  IAS  project  showed  that  it  could 
reliable  characterize  many  mine  explosions  and  underwater  explosions  based  on  the  observation  of 
spectral  modulations  due  to  ripple  fire.  We  have  observed  no  known  earthquakes  which  had 
spectral  modulations  as  strong  as  those  observed  from  known  ripple-fired  blasts,  detectable  by  the 
MHRSY  method,  which  suggests  that  the  false- alarm  rate  for  mine-blast  detection  would  be  low. 
However,  there  have  been  eases  where  MHRSY  failed  to  pick  up  a  weak  modulation,  caused  by  a 
small  time  delay.  There  have  also  been  known  blasts  which  did  not  produce  observable 
modulations,  either  because  they  were  not  ripple-fired  or  because  the  longest  time  delay  of  the 
npple  fire  never  exceeded  0.05  seconds,  the  shortest  delay  time 


showing  the  locations  of  events  which  we  presume  are  underwater  blasts. 


SPECTRA  FOR  EVENT  \  1  SPECTRA  FOR  EVENT  I  SPECTRA  FOR  EVENT 


FKiL  RK  21:  NORESS  spectra  and  cepstra  for  the  presumed  underwater  blasts  shown  on  the  map  in  Figure  20. 
Each  event  has  a  strong  consistent  modulation  in  the  0.4  to  0.6  second  quefrencv  range. 


We  have  decided  to  include  the  MF.RSY  process  in  the  ISF.IS  system  as  a  model-based 
reasoning  technique  for  identification  of  mine  blasts  and  underwater  blasts.  The  process  will 
effectively  identify  all  ripple  fired  events  w  ith  longest  delays  in  excess  of  0.05  seconds.  The  major 
features  of  the  ISHIS  version  of  MF.RSY  is  that  the  decision  rules  for  determining  if  the  cepstral 
peaks  are  significant  and  whether  they  are  at  the  same  quefrency  for  different  phases  will  be  coded 
in  the  expert  system  rather  than  in  the  actual  MFRSY  code  itself.  This  will  enable  us  to  change  the 
rules  in  order  to  to  make  them  more  sensitive  to  weak  modulations  without  unduly  increasing  the 
false  alarm  rate.  We  anticipate  that  the  ISHIS  version  of  MFRSY  will  identify  a  greater  number  of 
ripple-fired  blasts  than  the  IAS  version.  What  percentage  of  all  events  this  method  actually 
identifies  will  be  of  major  interest  in  the  evaluation  of  the  effectiveness  of  MF.RSY  and  the  ISEIS 
system  overall. 
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SECTION  4.0 


CASE-BASED  REASONING:  UTILITY  OF  DYNAMIC  TIME 

WARPING  (DTW)  FOR  REGIONAL  EVENT  CHARACTERIZATION 

4.1  INTRODUCTION 

In  this  section,  we  discuss  one  of  the  primary  algorithms  we  plan  to  implement  in  ISEIS 
for  case- based  reasoning.  We  have  observed  in  our  analysis  of  regional  coda  shapes,  in  the  form 
of  incoherent  beams  computed  for  different  frequency  bands,  that  we  can  often  identify  event  types 
from  the  distinctive  shape  of  the  coda.  Thus,  we  plan  to  implement  a  matching  algorithm  for 
characterizing  an  event  in  terms  of  the  similarity  of  its  coda  to  those  of  nearby  reference  events. 

Dynamic  time  warping  (DTW)  is  an  algorithm  for  matching  togv  cr  time-series  functions 
on  dissimilar  time  scales.  D  I'W  has  been  developed  in  the  speech  processing  field  for  word 
recognition.  In  this  application,  spoken  words  are  recognized  by  matching  a  signal  sequence  for 
the  unknow  n  word  against  a  set  of  template  patterns  for  previously  sampled  words.  Because  the 
unknown  word  probably  is  spoken  at  a  different  speed  than  any  of  the  template  words,  the  time 
scale  on  the  unknown  word  will  also  not  match  that  of  any  of  the  templates.  The  time  scales  of 
either  the  unknown  word  or  the  template  must  be  stretched  or  contracted  in  order  to  get  the  time 
scales  to  be  aligned  for  matching  purposes.  DTW  is  a  methodology  for  accomplishing  this 
rescaling  autom,-,’cally. 


The  idea  of  applying  DTW  to  seismic  signal  recognition  was  first  investigated  by  Anderson 
(  1<)M  )  Baumgardt  ( 19X7;  1990)  suggested  its  use  for  recognizing  regional  seismic  waveforms 
recorded  at  the  regional  seismic  array.  NORESS.  As  mentioned  above,  we  have  found  that 
regional  seismograms  recorded  at  NORESS  have  very  distinctive  shapes,  as  revealed  by 
incoherent  beam  waveform  envelopes.  It  has  been  possible,  for  example,  to  immediately  tell  the 
difference  between  western  Norway  mine  blasts  and  earthquakes  by  examining  the  character  of  the 
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waveform  shapes  of  the  events.  Baumgardt  (19X7;  1990)  suggested  matching  incoherent  beams 
using  1)TW  to  tjuan titat i vely  estimate  the  similarity  between  two  waveforms  for  discrimination 
purposes. 

Given  some  new  event  of  unknown  identity,  we  w  ish  to  find  all  the  previous  case  events 
which  are  most  similar  to  it,  in  terms  of  the  shape  of  its  waveform  in  a  set  of  prefiltered  frequency 
bands.  We  match  successively  the  incoherent  beams  or  templates  in  different  filter  bands  of  the 
unknow  n  event  against  the  templates  of  various  master  events  in  the  same  frequency  bands.  The 
figure-of-merit  for  the  matches,  which  is  the  distance  of  separation  in  log-rms  units  between  the 
unknown  event  and  the  reference  event,  provides  a  quantitative  measurement  of  how  similar  the 
unknown  event  is  to  the  reference  events.  In  this  application,  regional  events  that  need  to  be 
matched  may  not  occur  at  the  same  distance  from  the  array  as  the  reference  events.  Events  at 
different  distances  would  have  different  time  delays  between  the  various  phases  thus  necessitating 
an  algorithm  which  would  scale  the  unknown  event  and/or  the  template  event  to  a  common  time 
scale. 


This  section  is  a  brief  description  of  the  algorithm  and  our  initial  results  of  applying  the 
algorithm  on  a  set  of  earthquakes  and  explosions  in  western  Norway  recorded  at  the  NORESS 
array.  These  events  were  studied  during  a  research  project  relating  to  applying  case-based 
reasoning  to  seismic  event  identif  ication  (Baumgardt.  1990). 

4.2  DYNAMIC  TIME  WARPING  ALGORITHM 

Dynamic  time  warping  (DTW)  is  a  dynamic  programming  algorithm  used  to  find  the 
optimal  mapping  of  a  test  signal  time  axis  and  a  reference  signal  time  axis  onto  a  common  time 
axis.  The  approach  we  use  is  that  suggested  by  Anderson  ( 19X1 )  based  primarily  on  the  work  of 
Mvers  ( 19X0). 
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We  define  the  test  signal  as 


I  in  )  =  1 1  {n  1),/  ( n  2),  t  ( n  3),  ,.| 


and  the  reference  signal  is  defined  as 

Rin  )  =  \r  in  1  ),r  in  2),  r  in  3).  ...j 

where  t  and  n  represent  the  value  and  time  indices  of  the  test  signal  and  r  and  m  are  the  same  for 
the  reference  signal.  In  our  application,  /'and  R  would  be  the  incoherent  beams  in  some  filter 
frequency  band  for  the  new  event  to  be  identified  and  a  reference  event  of  known  ideniity.  For 
seismic  incoherent-beam  or  envelope  matching,  we  use  log  rms  values  computed  within  successive 
time  w  indow  s  in  the  seismogram  and  averaged  across  the  array.  Taking  the  logarithm  of  the  mis 
values  serves  to  enhance  some  of  the  envelope  shape  details  in  the  envelope  trace.  Also,  the 
difference  values  computed  by  DTW  are  in  log  units,  which  are  related  to  log  amplitude  values 
commonly  used  to  compute  seismic  magnitudes. 

A  warping  function  cij(k).  is  defined  which  matches  test  signal  index  i  with  reference 
signal  index  /.  or 

<  t  k  )  =  ■  n  ( /  1  .  in  ( /  ) ! 

‘■i  j 


where  k  is  the  new  time  axis.  Computing  this  function  can  be  thought  of  as  a  path  optimization 
problem,  as  illustrated  in  Figure  22.  The  time  axes  of  the  test  and  reference  series  represent  the 
horizontal  and  vertical  axes  of  the  grid  The  points  in  the  matrix  grid  represent  all  possible  match 
combinations  t/./i  for  the  time  points  of  the  two  time  functions.  It  we  just  lined  up  the  time  series 
and  matched  them  point  by  point  without  time  warping,  we  would  be  taking  a  diagonal  linear  path 
through  the  grid  in  Figure  22.  Warping  the  time  axes  involves  finding  a  non  linear,  non-diagonal 
path  through  the  grid  which  optimizes  some  match  function  between  the  two  time-series. 
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SCHEMATIC  ILLUSTRATION  OF  DTW  MATCH  AS  A  PATH 
MINIMIZATION  PROBLEM. 


FKiURE  22:  Schematic  illustration  of  DTW  match  as  a  path  minimization 
problem. 


The  optimal  path  through  the  grid  is  found  by  minimizing  the  accumulated  distance  function 


o 
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where  c(k)  is  the  warp  function  at  grid  point  k  (this  defines  the  i,j  values),  d  lc(k)l  is  the  local 
distance  measure  at  this  grid  point,  W(k)  is  a  local  weighting  function,  and  N(k)  is  the 
corresponding  normalization. 

Solving  the  DTW'  problem  requires  specification  of  the  following  parameters  (Myers,  et  al., 
1978):  (lj  Local  Distance  measure.  (2)  Axis  Orientation,  (3)  Endpoint  Constraints,  (4)  Local  Path 
Constraints,  and  (5)  Global  Path  Constraints.  In  addition,  W(k)  and  N.  the  weighting  function 
and  corresponding  normalization,  must  be  specified.  We  have  used  the  work  of  Anderson  ( 1981 ) 
and  Rabiner  and  Brow  n  ( 1982)  for  our  specifications  of  these  constraints. 


The  DTW  algorithm  used  for  seismic  template  matching  for  the  1SEIS  project  has  the 
following  specifications: 


Local  Distance  Measure 


The  local  distance  measure  for  seismic  template  matching  is  given  by  Anderson  (1981)  as 


d  Tic  Ik  )i  R  (  <  (  k  )  )=|  7  In  (  i  j)  K  (  rn  (  /  )) 


This  means  that  the  waveforms  are  aligned  based  on  the  area  under  their  normalized 
envelopes 
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Axis  Orientation 


Rabiner.  et  al  ( 1978)  report  best  pert'onnance  with  the  test  sequence  index  on  the  abscissa 
(independent  variable)  and  the  reference  sequence  index  on  the  ordinate  (dependent  variable),  as 
shown  in  Figure  22. 

Endpoint  Constraints 

We  require  the  path  to  obey  the  constraints  all  =1  and  c(M)  -  M.  This  assumes  the  end 
points  of  the  reference  and  test  signals  are  accurately  determined.  We  define  the  start  point  as  the 
first  phase  onset,  which  is  usually  the  Pn.  Pi;,  or  teleseismic  P  phase,  and  the  end  point  as  the 
place  w  here  the  Lf>  phase  drops  a  certain  percent  of  its  maximum.  This  assumption  simplifies  the 
optimum  path  problem. 

Local  Path  Constraints 

This  constraint  defines  the  movement  from  one  point  to  another.  Anderson  (1981 )  found 
that  the  type  II  d  matching  function  given  by  Myers  (1980)  works  reasonably  well  for  seismic 
envelope  matching.  For  this  algorithm  the  valid  paths  are: 

n  I i  ).  m  <  j )  ->  n  <i—l),  m  <j  -  I >. 


n  (i  ).  m  (./  )  ->  n  (i  -2),  m  (j--l  j. 
n(i).  m(j)  >  n<i—  I).  rn(j  2i. 

Figure  23  shows  these  paths.  This  is  a  modified  form  of  the  local  constraints  in  the 
symmetric  DTW  algorithm  for  a  maximum  slope  value  of  two  described  by  Itakura  ( 1*775 ).  This 
t>pe  of  movement  through  the  grid  allows  for  local  continuity  and  ensures  monotomcally 
increasing  time-series  indices  (Myers.  1980) 
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Global  Path  Constraints 

These  are  limitations  on  where  the  path  can  wander  in  the  grid.  The  slope  constraint  is  a 
natural  extension  of  the  endpoint  and  local  path  constraints.  An  additional  global  constraint  is 
added  which  keeps  the  warp  factor  between  the  two  envelopes  from  becoming  too  large.  It  is 
defined  by  R  where 

j  n  (  i  )  -  m  ( j  )  j  <  R 

The  combination  of  these  global  constraints  describe  a  parallelogram  on  the  grid  which 
encloses  all  the  points  which  could  possibly  occur  in  the  path,  as  shown  in  Figure  24. 

Weighting  Function  and  Corresponding  Normalization 

We  use  the  weighting  function  type  DTW  algorithm,  defined  by  Myers  (1980),  which  is 
W  i  k  )  =  n  (  i  )-«(/-  1  )  +  m  t  /  )  -  m  ( j  -  1 ) 

Myers  ( 19X0)  defines  a  normalization  function  which  is  independent  of  path  in  order  to 
solve  for  the  optimal  path  efficiently.  This  definition  is 

N  -  n  (i  )  +  m  (  /  ) 

Finding  the  minimum  for  the  accumulated  distance  function  can  be  done  by  calculating  the 
local  distance  at  every  grid  point  and  then  searching  all  possible  paths  from  (1,1)  to  (N,  M). 
Alternatively,  the  accumulated  distance  can  be  computed  using  a  dynamic  programming  technique 
(Myers,  1980).  The  dynamic  programming  principle  used  is  that  a  globally  optimal  path  is  locally 
optimal.  We  can  successively  build  the  optimal  path  from  smaller  optimal  paths  (Bellman  and 
Drevfus.  1962)  The  dynamic  programming  technique  and  the  constraints  reduces  significantly  the 
number  of  computational  operations  over  matching  by  complete  search  of  all  possible  paths.  Other 
techniques  arc  currently  being  researched  by  the  speech  processing  community  to  reduce  the 
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FKil  RF  24:  Slope  constraints  for  (he  paths  allowed  for  movement  through  the 
DTW  grid  shown  in  Figure  22. 


number  of  operations  for  the  grid  search  problem  (Brown  and  Rabiner,  1982;  Vidal  et  al,  1988). 
Other  research  has  been  directed  towards  implementing  the  DTW  algorithm  on  array  processors 
(e  g.,  Westc  et  al.  1983). 

The  DTW  algorithm  described  by  Myers  ( 1980)  computes  a  partial  accumulated  distance  at 
each  path  point  as  it  steps  through  the  grid.  The  partial  accumulated  distance  is  defined  recursively 

as 


l)  .  n  ( / 

1  ),  m  (  /  -  1  ) 
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where  is  the  local  distance  measure,  defined  above.  For  each  point  in  the  path  we  store  the  partial 
accumulated  distance,  l),\  \n  ( t  ).  m  (:  )  |.  and  the  index  of  the  previous  grid  point  in  the  path 
Then,  at  grid  point.  (\.  M)  the  total  accumulated  distance  is  given  by 

l)  =  Da  (N  ,  M  )/  N 

The  normalization  factor  does  not  affect  the  location  ot  the  optimal  path  since  it  is  defined 
as  path  independent  (Myers,  1980). 

The  output  accumulated  distance  value  is  a  measure  of  the  difference  between  the  test 
waveform  and  the  reference  waveform.  It  takes  into  account  the  time  alignment  of  the  waveforms 
based  on  the  area  under  their  normalized  envelopes.  This  DTW  "score"  can  now  be  used  in  a 
nearest  neighbor  decision  rule  to  find  the  best  matching  reference  template  or  it  can  be  used  to 
cluster  similar  events 

4.3  APPLICATION  TO  WKSTF.RN  NORWAY  EVENTS 

In  this  secnon.  we  discuss  the  application  of  this  algorithm  to  the  western  Norway  events 
discussed  m  detail  by  Baumgardt  (1990)  and  Baumgardt  and  Young  (1990).  The  events  are 
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shown  on  the  map  in  Figure  25  reproduced  from  Baumgardt  (1990).  The  events  we  originally 
believed  to  be  earthquakes  are  labeled  as  "Q"  on  the  map.  Flowever,  Baumgardt  and  Young 
(1990)  have  argued  that  many  of  the  "Q"  events  are  actually  blasts. 

The  event  we  chose  as  our  unknown  event  is  called  eqa  and  occurred  in  the  Stavanger 
region.  This  event  is  labeled  as  T.  '  on  the  map  in  Figure  25.  (Note:  This  event  is  referred  to  as 
event  L  in  Baumgardt  (1990)  and  was  changed  to  Q15  in  Baumgardt  and  Young  (1990).) 
Originally,  we  thought  that  this  event  was  an  earthquake  because  it  was  not  reported  as  being  a 
blast  in  the  Bergen  regional  bulletin.  However,  the  spectra  the  regional  phases  all  showed  the 
same  time  independent  spectra!  modulations,  which  appeared  to  be  produced  by  ripple-fire.  Also, 
our  studies  of  the  .S-to-P  amplitude  ratio  discriminant  revealed  that  the  event  has  ratios  more 
consistent  with  explosions  than  earthquakes.  Thus,  we  concluded  that  the  event  was  actually  an 
unannounced  blast. 

DTW  was  run  on  a  suite  of  incoherent  beams,  computed  for  eight  filter-frequency  bands, 
for  the  event  eqa  us  the  unknown  event  and  all  the  other  events  as  reference  events.  Figure  26 
shows  a  How  di.  -ram  'or  the  process  First,  the  incoherent  beams  are  lined  up  in  time.  In  this 
case,  we  lined  up  the  templates  on  the  first  arrival  Pn  phase.  Note  that  on  each  template,  a  phase 
identification  defines  a  time  interval,  from  the  first  break  of  the  phase  to  the  time  when  the  next 
phase  comes  in.  Then,  the  amplitudes  are  aligned  on  the  /.g  amplitude.  The  Lf>  phase  was  often, 
but  not  alwu>  s,  the  maximum  amplitude  phase  on  the  seismogram.  Finally,  all  the  points  after  Hn 
on  the  unknown  event  were  matched  up  to  reference  event  which  is  stretched  in  t.me  sufficient  to 
minimize  the  accumulated  distances,  as  explained  above. 

Figure  27  shows  a  DTW  solution  tor  the  match  of  the  incoherent  beam  of  this  event,  in  the 
X  to  16(1/  band,  compared  to  a  nearby  earthquake,  called  eq  9.  We  are  reasonably  sure  eq-9  was 
an  earthquake  because  the  local  news  media  in  Norway  reported  that  it  had  been  felt  (Svein 
Vlvkkeltveit  and  Frode  Ringdal.  personal  communication)  The  plot  in  Figure  27  shows  the  two 
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FKiliRK  25:  Map  showing  the  locations  of  the  events  used  in  this  study. 
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Mow  diagram  for  the  I)TV\  algorithm. 


74 


1.323 


A  c  c  d  i  s  t 
23.30 


A  v  e  d  i  s  t 


Pn 

P9 

S  n 
l9 


1. 22 
0.  H8 
0.37 
0.  OH 


Time  stretch  (from  PrO 

-3.00 
1.00 
6.  00 


FIGURE  27:  Example  of  a  DTW  solution  for  the  match  of  unknown  event,  eqa 
(top)  against  a  known  earthquake  eq^  (bottom). 
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incoherent  beams  for  the  events,  with  the  phase  identifications  indicated.  The  dashed  lines  indicate 
the  points,  which  mark  the  first  and  last  point  of  each  time  interval  assigned  to  each  phase,  that 
were  connected  by  the  DTW  algorithm.  At  the  bottom  of  the  plot  is  a  table  showing  the  total 
accumulated  distance,  the  average  distance  for  the  time  interval  assigned  to  each  phase,  and  the 
time  stretch,  which  is  how  much  the  time  of  each  phase  of  the  reference  event  (on  the  bottom)  had 
to  be  warped  for  the  envelopes  to  best  match.  To  the  right,  the  "Template  Relative  Magnitude" 
refers  to  the  apparent  difference  in  local  magnitude  based  on  the  difference  in  the  log-rms 
amplitudes  of  the  Lg  phases. 

This  result  shows  that  the  event  eqa  does  not  match  the  nearby  earthquake  very  well.  The 
biggest  mismatch  is  in  the  Pn  phase,  which  was  1.22  log  units  lower  for  eq-9  compared  with  our 
test  event,  eqa.  This  mismatch  suggests  that  event  eqa  produced  stronger  Pn  energy  relative  to  L\> 
than  the  earthquake.  —  .  -  -  ■  - 

Figure  28  shows  the  DTW  solution  which  gave  the  best  match  in  the  8  to  16  Hz  band,  an 
event  we  called  ql4.  Interestingly.  ql4  was  not  near  event  eqa  and  was  in  fact  316  km  to  the  north. 
However,  they  both  were  almost  exactly  the  same  distance  from  NORESS,  at  373  km.  Event  ql4 
occurred  due  west  of  NORESS  whereas  event  eqa  was  southwest.  The  DTW  accumulated 
distance  value  was  small  (8.33)  compared  to  the  earthquake  match  in  Figure  26.  Also,  visually, 
the  incoherent  beams  look  very  similar.  Baumgardt  and  Young  (1990)  also  concluded  that  event 
ql4  was  another  unannounced  blast,  like  eqa.  This  result  shows  that  the  presumed  unannounced 
blasts  have  very  similar  waveforms  and  that  their  similarity,  in  spite  of  the  fact  that  they  are  not  in 
the  same  region,  indicates  that  differences  in  propagation  path  do  not  have  a  great  effect  on  the 
waveform  characteristics  of  these  events. 


Figure  29  shows  the  comparison  of  the  event  with  a  known  mine  explosion,  ex2.  which 
occurred  at  the  Bl.A  site.  This  comparison  is  interesting  because  it  is  an  example  of  a  comparison 
of  two  events  at  very  different  distances:  BLA  is  at  301  km  whereas  eqa  was  at  373  km.  This  72 
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FIGURE  29:  DTW  solution  for  the  match  of  a  known  explosion,  ex  I,  to  the 
unknown  event,  eqa. 
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km  difference  in  distance  required  that  ex2  had  to  be  stretched  to  attain  the  best  match  to  eqa.  As 
shown  in  Figure  29,  the  Lg  of  ex2  had  to  be  warped  by  1 1  seconds  to  get  the  best  alignment. 
However,  after  doing  this,  the  accumulated  distance  was  8.53  log  rms  units,  which  was  the  second 
smallest.  Clearly,  after  time  warping,  eqa  more  closely  resembles  distant  mine  explosions  than 
nearby  earthquakes,  which  confirms  the  contention  of  Baumgardt  (1990)  and  Baumgardt  and 
Young  (1990),  based  on  other  evidence,  that  eqa  is  an  unannounced  blast. 

Figure  30  shows  plots  of  the  six  best  matching  events  to  eqa  in  the  2  to  4  Hz  band,  with  the 
incoherent  beams  plotted  up  in  order  of  the  accumulated  distance  of  separation.  Figure  31  shows 
the  same  display  for  the  incoherent  beams  in  the  8  to  16  Hz  band.  These  plots  show  that 
essentially  the  same  events  match  best  in  both  frequency  bands,  except  that  the  order  is  not  always 
the  same.  We  find  that,  for  all  the  frequency  bands,  the  events  which  match  best  were  either 
known  blasts  or  events  which  we  strongly  suspect  are  blasts,  based  on  other  evidence. 

4.4  CONCLUSIONS 

This  study  has  demonstrated  that  DTW  is  a  very  useful  method  for  quantitatively  finding 
regional  events  which  have  similar  waveforms.  DTW  has  been  implemented  in  1SEIS  for  the 
purpose  of  doing  case-based  reasoning  to  determine  how  similar  an  event  is  to  previous  events  in 
the  same  region.  Moreover,  we  also  plan  to  use  the  identity  of  the  best  matching  event  templates  to 
identify  the  new  event  if  the  cumulative  distance  of  separation  is  less  than  a  specified  maximum 
threshold.  For  example,  in  our  study,  if  we  set  the  minimum  separation  distance  at  10,  we  would 
find  from  Figure  31  that  the  best  matching  events  were  two  known  mine  blasts  (exl  and  ex2)  and 
one  strongly  suspected  blast  (ql4).  (Note:  The  event  also  matched  well  to  q!6,  with  a  distance  of 
10.7.  We  also  believe  that  ql6  is  an  unannounced  blast.)  We  also  note  that  all  known  earthquakes 
matched  with  distances  much  greater  than  10.  Thus,  we  would  conclude  that  event  eqa  is  most 
likely  a  blast  since  it  matches  the  explosion  templates  as  a  group  more  closely  than  the  earthquake 
templates. 
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FIGURE  30:  The  best  matching  templates  in  the  2  to  4  Hz  band  to  the  unknnwn 
event,  eqa,  shown  at  the  bottom.  a  a  to  tne  unknown 
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FIGURE  31:  The  best  matching  templates  in  the  8  to  16  Hz  band  to  the  unknown 
event,  eqa,  shown  at  the  bottom. 
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SECTION  5.0 


OVERALL  CONCLUSIONS 

This  report  has  discussed  some  of  the  design  and  computational  issues  related  to  the  1SE1S 
system.  In  this  section,  we  summarize  our  conclusions. 

System  Design 

( 1 )  The  ISEIS  approach  to  discrimination  will  be  to  test  each  individual  event  against  a 
matrix  of  discriminants,  rather  than  using  the  multivariate  decision  theory  approach 
to  discrimination.  This  approach  will  be  facilitated  by  the  spreadsheet  display, 
discussed  in  Section  1.0  (Figure  4),  where  each  column  of  the  spreadsheet  is  a 
single  discriminant  and  each  row  is  an  event  being  identified.  A  voting  scheme  will 
be  developed  where  each  discriminant  makes  an  individual  decision  about  the 
identity  of  the  event  at  some  confidence  level.  Then,  each  of  the  decisions  is 
combined  in  a  weighted  average  of  the  confidences  and  and  overall  decision  about 
the  events  identity  will  be  made. 

(2)  Discriminants  will  be  used  to  attempt  to  identify  the  event  as  a  specific  source  type, 
i.e.,  earthquake,  mine  blast,  underwater  explosion,  or  other  explosion.  Otherwise, 
the  event  will  be  classed  as  unidentified.  We  refer  to  this  type  of  event 
characterization  as  "model-based  reasoning,"  or  MBR.  We  recognize,  however, 
that  the  science  of  seismic  discrimination  has  not  yet  produced  a  completely  reliable 
method  for  identifying  seismic  events.  If  events  cannot  be  reliably  identified,  the 
events  can  still  be  compared  to  other  events  in  the  same  region  by  comparing  the 
waveform  features  of  the  event  to  be  identified  to  those  of  nearby  reference  events. 
We  call  this  kind  of  discrimination  "case-based  reasoning,"  or  CBR.  Moreover, 
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MBR  and  CBR  can  both  be  applied  in  a  complementary  manner  to  each  event  to 
completely  characterize  the  event. 


(3)  Certain  discriminants  can  be  used  both  in  MBR  and  CBR.  For  example,  the 
discriminants  discussed  in  Section  2.0  could  by  applied  in  both  ways.  Amplitude 
ratios  may  be  used  to  identify  events  as  a  specific  event  type.  Also,  they  can  be 
compared  with  amplitude  ratios  of  reference  events  in  CBR.  Other  discriminants, 
such  as  the  ripple-fire  detection  method  (MERSY)  discussed  in  Section  3.0, 
probably  only  work  as  only  as  MBR  methods.  However,  the  DTW  algorithm, 
described  in  Section  4.0,  although  primarily  a  waveform  pattern  matcher  and  thus  a 
CBR  method,  could  also  be  considered  an  MBR  technique  if  the  reference  events 
are  of  known  source  type.  Thus,  if  an  event  matched  best,  say,  to  mine-blast  coda 
shapes  in  the  same  region,  the  event  could  be  classed  as  a  mine  blast,  an  MBR  type 
classification,  even  though  a  CBR  matching  method  has  been  used  to  make  the 
decision. 


(4)  We  consider  effective  explanation  displays  to  be  an  important  feature  of  ISEIS. 
Users  must  be  able  to  get  adequate  explanations  of  the  processing  results  at  some 
level.  The  ISEIS  system  will  provide  different  levels  of  explanation,  where  the 
highest  level,  such  as  the  map  or  spreadsheet,  would  provide  broad,  summary 
views  of  the  event  characterization  results,  whereas  the  lowest  levels  would  give 
more  detailed  descriptions  (e.g.,  textual  descriptions,  graphs,  scatter  plots)  of  the 
processing  results. 

(5)  Because  of  the  incomplete  knowledge  of  seismic  discrimination  and  the  on-going 
research  in  this  field,  ISEIS  need  to  be  flexible  enough  to  be  easily  modified  to 
include  new  discriminants  and  to  change  or  remove  existing  discriminants.  To  this 
end,  ISEIS  is  being  designed  to  be  a  modular  system,  with  the  signal  analysis 
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processing  being  separated  from  the  discriminant  decision  processing.  This 
separation  will  be  effected  by  coding  the  decision  criteria  as  CLIPS  rules  in  a 
modifiable  ktjSwledge  base  of  an  expert  system.  Thus,  the  signal  analysis 
functions  would  compute  many  features  for  all  detected  phases  which  would  be 
stored  in  the  database,  but  the  rules  may  only  look  at  some  of  the  features.  In  fact, 
the  rules  themselves  will  be  treated  as  data  and  be  read  in  from  the  database  by  the 
system  along  with  waveform  feature  data.  After  more  data  has  been  analyzed, 
perhaps  using  some  of  the  low-level  displays  of  ISEIS,  other  features  than  those 
currently  reasoned  on  by  the  rules  may  be  found  that  are  more  effective.  We  may 
find  also  that  certain  discriminants  work  well  in  one  geographic  region,  but  less 
well  in  another,  perhaps  because  of  propagation  path  effects.  Moreover,  the 
threshold  values  for  discriminant  decisions  may  differ  from  one  region  to  another. 
Thus,  the  rules  can  be  modified  or  new  blocks  of  rules  can  be  added  in  the  form  of 
a  new  discriminant  or  premises  added  to  old  rules  to  make  them  region  specific. 
Because  the  rules  are  separated  from  the  signal  processing  algorithms  and  are 
treated  as  data,  changing  rules  or  adding  new  rules  will  not  require  a  complete 
change  and  recompilation  of  existing  code,  making  ISEIS  a  flexible  system  which 
can  easily  be  changed  as  new  discrimination  knowledge  is  obtained. 

Discriminants 

( 1 )  As  mentioned  above  and  in  Section  2.0,  our  knowledge  of  regional  discrimination 
is  at  best  problematic.  Thus  initially,  ISEIS  will  be  designed  to  include  a  set  of 
discriminants  which  look  promising  now  based  on  the  most  recent  research. 
However,  ISEIS  will  be  designed  to  be  flexible  enough  to  be  easily  modified  to 
take  advantage  of  the  results  of  new  research  in  the  field.  In  fact,  the  processing 
and  display  capabilities  of  ISEIS  and  its  access  to  expanded  databases  will  make  a 
useful  research  tool  in  seismic  discrimination. 
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(2)  The  initial  prototype  of  1SEIS  will  he  primarily  designed  to  process  small  events  at 
regional  seismic  distances.  Teleseismic  discrimination  will  be  less  emphasized, 
although  it  has  been  more  researched  and  may  be  important  in  event  identification  in 
third  countries  (i.e.,  not  the  Soviet  Union)  where  there  are  limited  in-country 
seismic  assets.  We  consider  it  beyond  the  scope  of  1SEIS  to  include  the  signal 
analysis  techniques  needed  for  teleseismic  discrimination,  such  as,  for  example 
computation  of  Ms  or  mb,  location  and  depth  estimation,  since  these  processes  will 
be  accomplished  by  the  front-end  NMRD  processes.  However,  ISEIS  will  have 
access  to  this  data  through  the  Oracle  database  and  can  include  rules  for  teleseismic 
discrimination.  Our  current  plan  is  to  include  rules  in  ISEIS  for  the  teleseismic 
discriminants,  Ms-mb  and  depth,  if  this  information  is  available  in  the  database. 

(3)  Based  on  our  research,  described  in  Section  2.0  and  those  of  other  workers  in  the 
field,  we  have  found  the  relative  excitation  of  regional  P  and  5  waves  to  provide  the 
greatest  discriminatory  power,  although  we  have  seen  regional  variations  in  these 
amplitudes  which  arc  undoubtedly  caused  by  lateral  variations  of  propagation 
effects.  Initially,  we  plan  to  include  at  least  two  discriminants  based  on  high- 
frequency  amplitude  ratios,  Pn/Sn  and  Pn/Lg.  The  initial  rules  will  be  designed  to 
best  identify  events  in  western  Norway  recorded  at  NORESS,  since  this  is  the 
region  for  which  we  currently  have  the  most  data.  As  ISEIS  is  run  on  more  data  in 
future,  these  discriminants  will  be  made  region  specific. 

(4)  Lg  spectral  ratio  proved  to  be  less  useful  as  a  discriminant,  for  currently  unknown 
reasons.  However,  ISEIS  will  still  be  designed  to  include  this  discriminant, 
although  it  may  be  weighted  less  than  other  discriminants  in  an  overall  decision 
about  the  identity  of  the  event.  Moreover,  spectral  ratio  computations  will  be  made 
on  all  other  regional  phases  and  will  be  available  for  display  for  research  purposes. 
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We  will  also  include  a  path-correction  scheme,  which  can  be  used  to  correct  spectral 
rados  for  the  effects  of  anelastic  attenuation  along  the  path.  Again,  a  regionalizadon 
scheme  is  being  developed  which  will  allow  region-specific  path  corrections  to  be 
included. 

(5)  We  believe  the  ripple-fire  detector  method,  MERSY,  to  be  a  very  useful  technique 
for  the  detection  of  non-nuclear,  chemical  blasts,  if  they  are  ripple-fired.  Moreover, 
blasting  underwater  may  be  also  identified,  based  on  spectral  modulations  produced 
by  bubble-pulse  interference.  Of  course,  blasts  underwater  could  still  be  nuclear, 
and  this  method  alone  cannot  insure  that  any  kind  of  blasting  is  not  nuclear.  The 
method  has  proven  to  be  effective  for  identifying  many  ripple-fired  chemical  and 
will  at  least  ensure  that  such  small  events  are  not  incorrectly  identified  as 
earthquakes  or  just  classed  as  unidentified. 

(61  DTW  will  be  implemented  as  our  primary  CBR  method.  However,  the  algorithm 
must  still  be  regarded  as  experimental,  although  our  initial  tests  of  the  technique 
have  been  promising.  The  method  can  only  be  applied  for  events  which  have 
nearby  reference  events  or  reference  events  at  comparable  distance.  Also,  DTW  can 
only  be  applied  to  events  which  have  relative  clear,  uncontaminated  codas.  Thus, 
DTW  could  not  actually  be  applied  to  events  which  have  interfering  events.  For 
this  reason,  adequate  displays  of  the  DTW  results  will  be  provided  to  the  user  to 
judge  the  validity  of  the  match  based  on  the  visual  appearance  of  the  codas 
themselves. 

(7)  For  contaminated  events,  other  case-based  methods  than  DTW  will  be  developed  to 
match  feature  measurement  to  those  of  reference  events.  Our  current  plan  is  to  use 
the  "script-matching'"  approach  described  by  Baumgardt  ( 1 DK7 )  and  Kandl  et  al 
(1987)  as  a  CBR  method  for  all  waveform  type  discriminants  (not  including 
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ME-.RSY ),  whore  we  compute  the  confidence  of  match  between  a  feature  for  the 
"current  event"  being  identified  and  average  feature  values  of  all  the  reference 
events  in  the  same  region.  A  regionalization  scheme  is  being  developed  in  which 
each  reference  event  will  be  tied  to  a  specific  reference  geographic  region.  The  new 
event  feature  will  then  be  matched  to  the  average  of  the  reference-event  features  for 
the  region  that  the  new  event  is  in  or  located  closest  to  it.  The  confidence  of  match 
will  then  be  reasoned  on  by  CBR  rules  to  determine  if  the  event  is  similar  or 
dissimilar  to  the  reference  events. 
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